/* ======================================================================== *\
!
! *
! * 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 18.05.2004  mazin@mppmu.mpg.de
// **********************************************************************************
// this macro is used to produce an alpha plot/plots for a chosen position in the sky map.
// In case of several subsamples, the overall alpha plot is also produced.
// The derotation is possible if Zd and Az are available for each event.
// Using of OFF data to estimate Nex and Significance is optional (if such data available)
//
// input: hillas parameter file/files
//      optional: a) ascii file with the estimated position of the source
//                b) ascii file with fit parameters from the OFF sample 
//                                  (needed for estimation of the background from OFF dada)
// output: alpha plot for each subsample, overall alpha plot
// **********************************************************************************


#define XSOUR   0.05  // [deg]
#define YSOUR   0.6  // [deg]
#define NUM     0


TString sourcename = "M87On";

const Bool_t ROTOPTION  = kFALSE;  // kFALSE: do not derotate, use camera coordinates
                                  // kTRUE: derotate into "quasi" sky coordinates

const Bool_t USEOFFDATA = kFALSE; // kFALSE: do not use OFF data to estimate the significance
                                  // kTRUE: use in addition OFF data to estimate the significance

const char * offfile = "paramOffCrab2701.dat"; // file with parameters of the OFF sample
                                               // needed only if USEOFFDATA = kTRUE

const Bool_t USEFILE = kFALSE; // kFALSE : use XSOUR, YSOUR 
                               // kTRUE: use "posfile" 
                               // to specify the position for which alpha plot has to be produced

const Bool_t SAMPLES = kFALSE; // kFALSE: one sample only
                               // kTRUE: several samples, usage of posfile is preferred!  

const char * posfile = "data/trackCrab1502_1305_Berlin_5deg.dat"; // file with best position
                                                                  // if USEFILE == kFALSE : not needed

TString dirname  = "~/data/M87/2004_02_17/";
TString filename = "M87HillasON.root"; //  NOTE in case of several subsamples (SAMPLES=kTRUE)
                                         // folllowing name construction is assumed:
                                         // file = dirname + sample + <number> + / + filename
                                         // e.g: ~/data/Crab/2004_01_27/sample2/CrabHillasON.root
#define tolerance 1e-3
/* ******************************************************   */
      /*  dynamical cuts  Crab 27th Jan 2004 */
/*
#define LENGTHMINParA 0.136  // deg
#define LENGTHMINParB 0.036  //
#define LENGTHMINParC -0.0038  //
#define LENGTHMAXParA 0.332  // deg
#define LENGTHMAXParB 0.037  //
#define LENGTHMAXParC 0.0261 //
#define WIDTHMINParA  0.063  // deg
#define WIDTHMINParB  0.013  //
#define WIDTHMINParC  0.0003  //
#define WIDTHMAXParA  0.123  // deg
#define WIDTHMAXParB  0.019  //
#define WIDTHMAXParC  0.0005 //
#define DISTMINParA   0.6   // deg
#define DISTMINParB   0.059  //
#define DISTMINParC   0.  //
#define DISTMAXParA   1.25   // deg
#define DISTMAXParB   0.059  //
#define DISTMAXParC   0.  //
*/

      /*  dynamical cuts  Mrk 421 and Crab 15th Feb 2004 */

#define LENGTHMINParA 0.12  // deg
#define LENGTHMINParB 0.034  //
#define LENGTHMINParC 0.  //
#define LENGTHMAXParA 0.32  // deg
#define LENGTHMAXParB 0.034  //
#define LENGTHMAXParC 0. //
#define WIDTHMINParA  0.055  // deg
#define WIDTHMINParB  0.013  //
#define WIDTHMINParC  0.0  //
#define WIDTHMAXParA  0.12  // deg
#define WIDTHMAXParB  0.013  //
#define WIDTHMAXParC  0.0 //
#define DISTMINParA   0.6   // deg
#define DISTMINParB   0.059  //
#define DISTMINParC   0.  //
#define DISTMAXParA   1.25   // deg
#define DISTMAXParB   0.059  //
#define DISTMAXParC   0.  //


#define SIZEMIN 2000
#define LEAKMAX 0.25
#define ASYMMIN  -0.1 

#include "mtools.C"

#define ALPHAMAX 15.

#define ALOFFMAX 90.
#define ALOFFMIN 30.

// histogram to store the sum of alpha plots
   TH1F histalphaAll("alpha plot", "alpha plot, size cut 2000, 15 deg", 36, 0., 90.);
   histalphaAll.SetXTitle("alpha (deg)");
   histalphaAll.SetYTitle("Counts");
   histalphaAll.SetDirectory(NULL);
   histalphaAll.SetFillStyle(4000);
   histalphaAll.UseCurrentStyle();


void alphaplots(Int_t num, Double_t XSOURCE, Double_t YSOURCE)
{

// for USEOFFDATA = kTRUE the ascii file with fit parameters of the OFF sample is read.

   if (USEOFFDATA == kTRUE)
   {
	    Int_t event, dummy, numbinsoff;
	    Float_t binwidthOff, xpos, ypos, aparOff, bparOff, chi2par;

	    FILE *fp;
	    fp = fopen(offfile,"r");    
	    while(fscanf(fp,"%d %f %f %f %f %f %f %f %f %f",
		       &event, &xpos, &ypos,  &aparOff, &bparOff, &chi2par,
		       &binwidthOff, &dummy, &dummy, &dummy) != EOF)
	    {
//cout << event << " " << xpos << " " << ypos << endl;
		     if(TMath::Abs(xpos-XSOURCE) < tolerance && TMath::Abs(ypos-YSOURCE) < tolerance) break;
	    }
	    fclose(fp);
   }


    gStyle->SetCanvasBorderMode(0);
    gStyle->SetCanvasBorderSize(0);
    gStyle->SetCanvasColor(10);
    gStyle->SetPadBorderMode(0);
    gStyle->SetPadBorderSize(0);
    gStyle->SetPadColor(10);
    gStyle->SetOptFit(1);
    gStyle->SetStatColor(10);
//    gStyle->SetOptStat(0);
    gStyle->SetPalette(1,0);



// the name of the Hillas parameter file, which has to be read in
    TString ddumy;

    if (SAMPLES == kFALSE)
        filename = dirname + filename;
    else
    {
       ddumy = dirname;
       ddumy += "sample"; 
       ddumy += num; 
       filename = ddumy + filename; 
    }

cout << "file to read :" << filename << endl;


// create histograms   
    MBinning bins;

   TH1F histlength;
   histlength.SetName("Length");
   histlength.SetTitle("Length");
   histlength.SetXTitle("Length [deg]");
   histlength.SetYTitle("Counts");
   histlength.SetDirectory(NULL);
   histlength.SetFillStyle(4000);
   histlength.UseCurrentStyle(); 


    bins.SetEdges(100, 0., 1.);
    bins.Apply(histlength);


   TH1F histwidth;
   histwidth.SetName("Width");
   histwidth.SetTitle("Width");
   histwidth.SetXTitle("Width [deg]");
   histwidth.SetYTitle("Counts");
   histwidth.SetDirectory(NULL);
   histwidth.SetFillStyle(4000);
   histwidth.UseCurrentStyle(); 

    bins.SetEdges(100, 0., 0.5);
    bins.Apply(histwidth);


   TH1F histsize;
   histsize.SetName("Size");
   histsize.SetTitle("Size");
   histsize.SetXTitle("Size");
   histsize.SetYTitle("Counts");
   histsize.SetDirectory(NULL);
   histsize.SetFillStyle(4000);
   histsize.UseCurrentStyle(); 

//    bins.SetEdges(100, 100., 2e5);
//    bins.Apply(histsize);

    bins.SetEdgesLog(100, 100., 10e5);
    bins.Apply(histsize);

   TH1F histalpha;
   histalpha.SetName("Alpha");
   histalpha.SetTitle("Alpha");
   histalpha.SetXTitle("Alpha [deg]");
   histalpha.SetYTitle("Counts");
   histalpha.SetDirectory(NULL);
   histalpha.SetFillStyle(4000);
   histalpha.UseCurrentStyle(); 

    bins.SetEdges(100, -100., 100.);
    bins.Apply(histalpha);

   TH1F histdist;
   histdist.SetName("Dist");
   histdist.SetTitle("Dist");
   histdist.SetXTitle("Dist [deg]");
   histdist.SetYTitle("Counts");
   histdist.SetDirectory(NULL);
   histdist.SetFillStyle(4000);
   histdist.UseCurrentStyle(); 

    bins.SetEdges(100, 0., 2.);
    bins.Apply(histdist);

   TH1F histmeanx;
   histmeanx.SetName("MeanX");
   histmeanx.SetTitle("MeanX");
   histmeanx.SetXTitle("MeanX [deg]");
   histmeanx.SetYTitle("Counts");
   histmeanx.SetDirectory(NULL);
   histmeanx.SetFillStyle(4000);
   histmeanx.UseCurrentStyle(); 

    bins.SetEdges(100, -1.8, 1.8);
    bins.Apply(histmeanx);

   TH1F histmeany;
   histmeany.SetName("MeanY");
   histmeany.SetTitle("MeanY");
   histmeany.SetXTitle("MeanY [deg]");
   histmeany.SetYTitle("Counts");
   histmeany.SetDirectory(NULL);
   histmeany.SetFillStyle(4000);
   histmeany.UseCurrentStyle(); 


    bins.SetEdges(100, -1.8, 1.8);
    bins.Apply(histmeany);

   TH1F histalphafinal;
   histalphafinal.SetName("ALPHA");
   histalphafinal.SetTitle("ALPHA");
   histalphafinal.SetXTitle("alpha [deg]");
   histalphafinal.SetYTitle("Counts");
   histalphafinal.SetDirectory(NULL);
   histalphafinal.SetFillStyle(4000);
   histalphafinal.UseCurrentStyle(); 

    bins.SetEdges(36, 0.0, 90.);
    bins.Apply(histalphafinal);

   TH1F histAssym;
   histAssym.SetName("Assymetry");
   histAssym.SetTitle("Assymetry");
   histAssym.SetXTitle("Assymetry");
   histAssym.SetYTitle("Counts");
   histAssym.SetDirectory(NULL);
   histAssym.SetFillStyle(4000);
   histAssym.UseCurrentStyle(); 

    bins.SetEdges(100, -1, 1.);
    bins.Apply(histAssym);

   TH1F histAssymM3;
   histAssymM3.SetName("Assymetry 3M");
   histAssymM3.SetTitle("Assymetry 3rd moment");
   histAssymM3.SetXTitle("Assymetry 3rd moment");
   histAssymM3.SetYTitle("Counts");
   histAssymM3.SetDirectory(NULL);
   histAssymM3.SetFillStyle(4000);
   histAssymM3.UseCurrentStyle(); 

    bins.SetEdges(100, -1., 1.);
    bins.Apply(histAssymM3);

   TH2F hist2xy("CoG","Center of Gravity", 100, -1.8, 1.8, 100, -1.8, 1.8);
   hist2xy.SetXTitle("MeanX [deg]");
   hist2xy.SetYTitle("MeanY [deg]");
   hist2xy.SetDirectory(NULL);
   hist2xy.SetFillStyle(4000);
   hist2xy.UseCurrentStyle();
 
   TH1F histLoverS;
   histLoverS.SetName("LoverS");
   histLoverS.SetTitle("LoverS");
   histLoverS.SetXTitle("LoverS");
   histLoverS.SetYTitle("Counts");
   histLoverS.SetDirectory(NULL);
   histLoverS.SetFillStyle(4000);
   histLoverS.UseCurrentStyle(); 

    bins.SetEdges(100, -0., 0.0006);
    bins.Apply(histLoverS);

   TH2F hist2lw("Length-Width", "correlation Length-Width", 100, 0.0, 1.0, 100, 0.0, 0.5);
   hist2lw.SetXTitle("Length [deg]");
   hist2lw.SetYTitle("Width [deg]");
   hist2lw.SetDirectory(NULL);
   hist2lw.SetFillStyle(4000);
   hist2lw.UseCurrentStyle();
 
   TH2F hist2lalpha("Length-Alpha", "correlation Length-Alpha", 100, 0.0, 1.0, 100, -100., 100.);
   hist2lalpha.SetXTitle("Length [deg]");
   hist2lalpha.SetYTitle("Alpha [deg]");
   hist2lalpha.SetDirectory(NULL);
   hist2lalpha.SetFillStyle(4000);
   hist2lalpha.UseCurrentStyle();
 
   TH2F hist2ldist("Length-Dist","correlation Length-Dist", 100, 0.0, 1.0, 100, 0.0, 1.7);
   hist2ldist.SetXTitle("Length [deg]");
   hist2ldist.SetYTitle("Dist [deg]");
   hist2ldist.SetDirectory(NULL);
   hist2ldist.SetFillStyle(4000);
   hist2ldist.UseCurrentStyle();
 
   TH2F hist2walpha("Width-Alpha","correlation Width-Alpha", 100, 0.0, 0.5, 100, -100., 100.);
   hist2walpha.SetXTitle("Width [deg]");
   hist2walpha.SetYTitle("Alpha [deg]");
   hist2walpha.SetDirectory(NULL);
   hist2walpha.SetFillStyle(4000);
   hist2walpha.UseCurrentStyle();
 
   TH2F hist2wdist("Width-Dist","correlation Width-Dist", 100, 0.0, 0.5, 100, 0.0, 1.7);
   hist2wdist.SetXTitle("Width [deg]");
   hist2wdist.SetYTitle("Dist [deg]");
   hist2wdist.SetDirectory(NULL);
   hist2wdist.SetFillStyle(4000);
   hist2wdist.UseCurrentStyle();
 
   TH2F hist2alphadist("Alpha-Dist","correlation Alpha-Dist", 100, -100., 100, 100, 0.0, 1.7);
   hist2alphadist.SetXTitle("Alpha [deg]");
   hist2alphadist.SetYTitle("Dist [deg]");
   hist2alphadist.SetDirectory(NULL);
   hist2alphadist.SetFillStyle(4000);
   hist2alphadist.UseCurrentStyle();
 
   TH1F histphi;
   histphi.SetName("TelPhia");
   histphi.SetTitle("Telescope Phi");
   histphi.SetXTitle("Phi [rad]");
   histphi.SetYTitle("Counts");
   histphi.SetDirectory(NULL);
   histphi.SetFillStyle(4000);
   histphi.UseCurrentStyle(); 

    bins.SetEdges(100, -10, 10);
    bins.Apply(histphi);

   TH1F histtheta;
   histtheta.SetName("TelTheta");
   histtheta.SetTitle("Telescope Theta");
   histtheta.SetXTitle("Theta [rad]");
   histtheta.SetYTitle("Counts");
   histtheta.SetDirectory(NULL);
   histtheta.SetFillStyle(4000);
   histtheta.UseCurrentStyle(); 

    bins.SetEdges(100, -2, 2.);
    bins.Apply(histtheta);

   TH1F aftercuthistlength;
   aftercuthistlength.SetName("Length");
   aftercuthistlength.SetTitle("Length");
   aftercuthistlength.SetXTitle("Length [deg]");
   aftercuthistlength.SetYTitle("Counts");
   aftercuthistlength.SetDirectory(NULL);
   aftercuthistlength.SetFillStyle(4000);
   aftercuthistlength.UseCurrentStyle();


    bins.SetEdges(100, 0., 1.);
    bins.Apply(aftercuthistlength);

   TH1F aftercuthistwidth;
   aftercuthistwidth.SetName("Width");
   aftercuthistwidth.SetTitle("Width");
   aftercuthistwidth.SetXTitle("Width [deg]");
   aftercuthistwidth.SetYTitle("Counts");
   aftercuthistwidth.SetDirectory(NULL);
   aftercuthistwidth.SetFillStyle(4000);
   aftercuthistwidth.UseCurrentStyle();

    bins.SetEdges(100, 0., 0.5);
    bins.Apply(aftercuthistwidth);

   TH1F aftercuthistsize;
   aftercuthistsize.SetName("Size");
   aftercuthistsize.SetTitle("Size");
   aftercuthistsize.SetXTitle("Size [photons]");
   aftercuthistsize.SetYTitle("Counts");
   aftercuthistsize.SetDirectory(NULL);
   aftercuthistsize.SetFillStyle(4000);
   aftercuthistsize.UseCurrentStyle();

    bins.SetEdgesLog(100, 100., 10e5);
    bins.Apply(aftercuthistsize);

   TH1F aftercuthistalpha;
   aftercuthistalpha.SetName("Alpha");
   aftercuthistalpha.SetTitle("Alpha");
   aftercuthistalpha.SetXTitle("Alpha [deg]");
   aftercuthistalpha.SetYTitle("Counts");
   aftercuthistalpha.SetDirectory(NULL);
   aftercuthistalpha.SetFillStyle(4000);
   aftercuthistalpha.UseCurrentStyle();

    bins.SetEdges(20, 0., 100.);
    bins.Apply(aftercuthistalpha);

   TH1F aftercuthistdist;
   aftercuthistdist.SetName("Dist");
   aftercuthistdist.SetTitle("Dist");
   aftercuthistdist.SetXTitle("Dist [deg]");
   aftercuthistdist.SetYTitle("Counts");
   aftercuthistdist.SetDirectory(NULL);
   aftercuthistdist.SetFillStyle(4000);
   aftercuthistdist.UseCurrentStyle();

   bins.SetEdges(100, 0., 2.);
    bins.Apply(aftercuthistdist);

   TH1F aftercuthistmeanx;
   aftercuthistmeanx.SetName("MeanX");
   aftercuthistmeanx.SetTitle("MeanX");
   aftercuthistmeanx.SetXTitle("MeanX [deg]");
   aftercuthistmeanx.SetYTitle("Counts");
   aftercuthistmeanx.SetDirectory(NULL);
   aftercuthistmeanx.SetFillStyle(4000);
   aftercuthistmeanx.UseCurrentStyle();

    bins.SetEdges(100, -1.8, 1.8);
    bins.Apply(aftercuthistmeanx);

   TH1F aftercuthistmeany;
   aftercuthistmeany.SetName("MeanY");
   aftercuthistmeany.SetTitle("MeanY");
   aftercuthistmeany.SetXTitle("MeanY [deg]");
   aftercuthistmeany.SetYTitle("Counts");
   aftercuthistmeany.SetDirectory(NULL);
   aftercuthistmeany.SetFillStyle(4000);
   aftercuthistmeany.UseCurrentStyle();


    bins.SetEdges(100, -1.8, 1.8);
    bins.Apply(aftercuthistmeany);

   TH2F aftercuthist2xy("CoG","Center of Gravity", 100, -1.8, 1.8, 100, -1.8, 1.8);
   aftercuthist2xy.SetXTitle("MeanX [deg]");
   aftercuthist2xy.SetYTitle("MeanY [deg]");
   aftercuthist2xy.SetDirectory(NULL);
   aftercuthist2xy.SetFillStyle(4000);
   aftercuthist2xy.UseCurrentStyle();



   const Int_t n = 100;
   Double_t binsize[n];

   Float_t nmin = 100.;
   Float_t nmax = 1e7;

   for(Int_t i=0; i<n; i++)
   {
    binsize[i] = pow(10., log10(nmin) + i * (log10(nmax) - log10(nmin)) / (n-1.));
   }

   TH2F hist2wsize("Width-Size", "correlation Width-Size", 100, 0.0, 0.5, n-1, binsize);
   hist2wsize.SetXTitle("Width [deg]");
   hist2wsize.SetYTitle("Size");
   hist2wsize.SetDirectory(NULL);
   hist2wsize.SetFillStyle(4000);
   hist2wsize.UseCurrentStyle();


 
   TH2F hist2alphasize("Alpha-Size","correlation Alpha-Size", 100, -100., 100., n-1, binsize);
   hist2alphasize.SetXTitle("Alpha [deg]");
   hist2alphasize.SetYTitle("Size");
   hist2alphasize.SetDirectory(NULL);
   hist2alphasize.SetFillStyle(4000);
   hist2alphasize.UseCurrentStyle();
 
   TH2F hist2distsize("Dist-Size","correlation Dist-Size", 100, 0.0, 1.7, n-1, binsize);
   hist2distsize.SetXTitle("Dist [deg]");
   hist2distsize.SetYTitle("Size");
   hist2distsize.SetDirectory(NULL);
   hist2distsize.SetFillStyle(4000);
   hist2distsize.UseCurrentStyle();
// end create histograms


 
    //
    // Now setup the tasks and tasklist:
    // ---------------------------------
    //

    MParList  plist;

    MTaskList tlist;
    plist.AddToList(&tlist);


    MReadMarsFile read("Events");
    read.DisableAutoScheme();

    read.AddFile(filename);

    MHillas mhillas;
    plist.AddToList(&mhillas);

    MHillasSrc mhillassrc;
    plist.AddToList(&mhillassrc);

    MHillasExt mhillasext;
    plist.AddToList(&mhillasext);

    MNewImagePar mnewimpar;
    plist.AddToList(&mnewimpar);

    MGeomCamMagic cam;
    plist.AddToList(&cam);
 
    MMcEvt mcevt;
    plist.AddToList(&mcevt);
 
    MPointingPos mpoint;
    plist.AddToList(&mpoint);

    MObservatory observ;
    plist.AddToList(&observ);

//    MRawRunHeader header;
//    plist.AddToList(&header);
    
    tlist.AddToList(&read);
   
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);

   if (!tlist.PreProcess(&plist))
      return;

    Float_t fMm2Deg = cam->GetConvMm2Deg();
    Int_t event = 0;
    Int_t filenumber = 0;

    Float_t ftheta, fphi, flength, fwidth, fsize, fmeanx, fmeany, falpha, fdist;
    Float_t fsingam, fcosgam;
    Double_t  xsournew, ysournew;
    Float_t fdelta, fleak, fconc1, fcosda, fassym, fassymM3;
    Int_t AsGrNull=0, AsLessNull=0;
    Int_t AsGrNullAfter=0, AsLessNullAfter=0;
    Float_t logsize, lgsize, lgsize2, tanbeta, beta;
    const Float_t LOG3000 = log(3000.);
    Char_t stringtriv1[80], stringlima[80], stringNex[80], stringsig[80];
    Char_t stringNexOnOff[80], stringLiMaOnOff[80];

// initial values:
     Float_t  xsource = XSOURCE;
     Float_t  ysource = YSOURCE;


    while (tlist.Process())
    {
    event++;

    if (mhillas->GetLength() != -1.)
     {
      // parameters:
      flength = (mhillas->GetLength()) * fMm2Deg;
      fwidth = (mhillas->GetWidth())*fMm2Deg;
      fsize = mhillas->GetSize();
      fmeanx = (mhillas->GetMeanX())*fMm2Deg;
      fmeany = (mhillas->GetMeanY())*fMm2Deg;
      falpha = mhillassrc->GetAlpha(); 
      fdist = (mhillassrc->GetDist())*fMm2Deg;
      fdelta = mhillas->GetDelta();
      fconc1 = (mnewimpar->GetConc1());
      fleak = mnewimpar->GetLeakage1();


//      ftheta = mcevt->GetTelescopeTheta();
      ftheta = mpoint->GetZd();
//      fphi = mcevt->GetTelescopePhi();
      fphi = mpoint->GetAz();
// cout << " phi : " << fphi << " theta : " << ftheta << endl;
      observ.RotationAngle(ftheta, fphi, fsingam, fcosgam);

      fassym = (mhillasext->GetAsym()) * fMm2Deg;
      fassymM3 = (mhillasext->GetM3Long()) * fMm2Deg;
      fcosda = mhillassrc->GetCosDeltaAlpha();

      if ((fassymM3*TMath::Sign(1.,fcosda)) > 0.) AsGrNull++;
      else AsLessNull++;

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

      // basic plots:

//    if (fsize > 3000.)
    if (fsize > 0.)
     {
      histphi.Fill(fphi,1.);
      histtheta.Fill(ftheta,1.);

      histlength.Fill(flength,1.);
      histwidth.Fill(fwidth,1.);
      histsize.Fill(fsize,1.);
      histLoverS.Fill(flength/fsize,1.);
      histmeanx.Fill(fmeanx,1.);
      histmeany.Fill(fmeany,1.);
      histalpha.Fill(falpha,1.);
      histdist.Fill(fdist,1.);
      hist2xy.Fill(fmeanx, fmeany, 1.);
     }

// some cuts:
     if (flength > 0.1 && flength < 0.32)
        if (fwidth > 0.06 && fwidth < 0.15)
          if (fdist > 0.6 && fdist < 1.3)
            if (fsize > 3000.)
//              if(sqrt(fmeanx*fmeanx + fmeany*fmeany) < 1.1)  
//              if((fassymM3*fcosda < 0.3 && fassymM3*fcosda > 0.02) || 
//                 (fassymM3*fcosda < -0.02 && fassymM3*fcosda > -0.2) )
//                  if(fassym*fcosda > 0.) 
                      {
                       histAssymM3.Fill(fassymM3*TMath::Sign(1.,fcosda), 1.);
                       histAssym.Fill(fassym*TMath::Sign(1.,fcosda), 1.);
                       if ((fassymM3*TMath::Sign(1.,fcosda)) > 0.) AsGrNullAfter++;
                       else AsLessNullAfter++;
                      }

// **********************************************************************  //
// calculate alpha and dist according to the source location:
           tanbeta = (fmeany - ysournew) / (fmeanx - xsournew);
           beta = TMath::ATan(tanbeta);
           falpha = (fdelta - beta) * 180./ TMath::Pi();
           fdist = sqrt((fmeany - ysournew) * (fmeany - ysournew) +
                             (fmeanx - xsournew) * (fmeanx - xsournew));

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

// **********************************************************************  //



      if (fsize > 3000.)
      {
      // correlations:
      hist2lw.Fill(flength, fwidth, 1.);
      hist2lalpha.Fill(flength, falpha, 1.);
      hist2ldist.Fill(flength, fdist, 1.);
      hist2walpha.Fill(fwidth, falpha, 1.);
      hist2wdist.Fill(fwidth, fdist, 1.);
      hist2alphadist.Fill(falpha, fdist, 1.);
      hist2wsize.Fill(fwidth, fsize, 1.);
      hist2alphasize.Fill(falpha, fsize, 1.);
      hist2distsize.Fill(fdist, fsize, 1.);
      }     

     // cuts:
//cout << " before the cuts"  << "size :" << fsize << endl;
                 logsize = log(fsize);
                 lgsize = logsize-LOG3000;
                 lgsize2 = lgsize*lgsize;
                 if ( fsize > SIZEMIN )
                 if ( fleak < LEAKMAX )
                 if ( flength > (LENGTHMINParA + LENGTHMINParB*lgsize + LENGTHMINParC*lgsize2) &&
                      flength < (LENGTHMAXParA + LENGTHMAXParB*lgsize + LENGTHMAXParC*lgsize2))
                 if ( fwidth > (WIDTHMINParA + WIDTHMINParB*lgsize + WIDTHMINParC*lgsize2)  &&
                      fwidth < (WIDTHMAXParA + WIDTHMAXParB*lgsize + WIDTHMAXParC*lgsize2) )
                 if ( fdist > (DISTMINParA + DISTMINParB*lgsize + DISTMINParC*lgsize2)  &&
                      fdist < (DISTMAXParA + DISTMAXParB*lgsize + DISTMAXParC*lgsize2) )
//                 if ((fassym*TMath::Sign(1.,fcosda)) > ASYMMIN)   // asymmcut
                 {
                         falpha = TMath::Abs(falpha);
                         histalphafinal.Fill(falpha,1.);
                         histalphaAll.Fill(falpha,1.);

                         aftercuthistlength.Fill(flength,1.);
                         aftercuthistwidth.Fill(fwidth,1.);
                         aftercuthistsize.Fill(fsize,1.);
                         aftercuthistmeanx.Fill(fmeanx,1.);
                         aftercuthistmeany.Fill(fmeany,1.);
                         aftercuthistalpha.Fill(falpha,1.);
                         aftercuthistdist.Fill(fdist,1.);
                         aftercuthist2xy.Fill(fmeanx, fmeany, 1.);

                 }


     }
    else filenumber++;
    }

//    cout << " conversion factor is: " << fMm2Deg <<  endl;
    cout << " events read in from file : " << event <<  endl;
    cout << " runs found in the file : " << filenumber <<  endl;

    Int_t startbinoff;
    Float_t Nex, Non, Noff, Sign, SignLiMa;
    Float_t normf, integon, integoff, NexOnOff, NoffOFF, SignOnOff, SignLiMaOnOff;
    Float_t binwidth = histalphafinal.GetBinWidth(1);
    Float_t numbinMax = ALPHAMAX/binwidth;
 
// **********************************************************************  //
/* fit parabel from 30 to 90 degrees   */
    TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", ALOFFMIN, ALOFFMAX);
    fitbgpar->SetLineColor(2);

    histalphafinal.Fit("fbgpar","WNR");

    Double_t apar = fitbgpar->GetParameter(0);
    Double_t bpar = fitbgpar->GetParameter(1);

    TF1 * bgoff = new TF1("bgoffON", parabfunc, 0., 90., 3);
    bgoff->SetParameters(apar, bpar, 1.);
    bgoff->FixParameter(0, apar);
    bgoff->FixParameter(1, bpar);
    bgoff->FixParameter(2, 1.);
    bgoff->SetLineColor(9);

/* end of the fit parabel from 30 to 90 degrees*/
// **********************************************************************  //



    if (!tlist.PostProcess())
      return;

    
    gStyle->SetOptStat(11);

// calculate significance: DO NOT USE FIT FOR Non!!!
    Non = 0.;
    for(Int_t i=1; i<=numbinMax;i++) Non += histalphafinal.GetBinContent(i);

    Noff = (1./3. * (fitbgpar->GetParameter(0)) * pow(ALPHAMAX,3.) +
           (fitbgpar->GetParameter(1)) * ALPHAMAX) /  binwidth;
    Nex = Non - Noff;

    Sign = Nex / sqrt(Nex + 2.* Noff);
 
    cout << " Non : " << Non << " Noff : " << Noff << " Nex : " << Nex << endl;
    cout << " significance : " << Sign << " sigma" << endl;

    SignLiMa = LiMa17(Non,Noff,1.);
    cout << " significance Li and Ma (17): " << SignLiMa << " sigma" << endl;

    Char_t stringsig[80];
 
    sprintf(stringsig,"S = %.2f sigma", Sign);
    sprintf(stringtriv1,"Signif:  S = %.2f sigma", Sign);
    sprintf(stringlima,"Li&Ma 17: S = %.2f sigma", SignLiMa);
    sprintf(stringNex,"N excess: Nex = %.d ", Nex);

// **********************************************************************  //
//  use OFF data to estimate background   *******************************  //
    TF1 * bgoff2 = new TF1("bgoffOFF", parabfunc, 0., 90., 3);
  if (USEOFFDATA == kTRUE)
  {
        // ON:

           integon = 0.;  // number of events between 30 and 90 degrees
	   numbinsoff = TMath::Nint((ALOFFMAX - ALOFFMIN)/binwidth);
           startbinoff = TMath::Nint(ALOFFMIN/binwidth) + 1;

           for (Int_t ik = 0; ik < numbinsoff; ik++)
           {
            integon += histalphafinal.GetBinContent(startbinoff+ik);  
           }
           // OFF:

           integoff = ((1./3. * aparOff * pow(90.,3.) + bparOff * 90.) -
                      (1./3. * aparOff * pow(30.,3.) + bparOff * 30.)) / binwidthOff;

           normf = integoff / integon;

           NoffOFF = (1./3. * aparOff * pow(ALPHAMAX,3.) +
                   (bparOff * ALPHAMAX)) /  binwidthOff / normf;

           NexOnOff = Non - NoffOFF;

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

        // calculate according to Li Ma:
           SignLiMaOnOff = LiMa17(Non,NoffOFF*normf,1./normf);

cout << " integon: " << integon << ",  integoff : " <<  integoff << ", normf: " << normf
     << ", NoffOFF : " << NoffOFF << ", Non : " << Non
     << ", NexOnOff : " << NexOnOff << ", SignOnOff : " << SignOnOff << endl;
 cout << " significance (LiMa 17): " << SignLiMaOnOff << " sigma" << endl;

 sprintf(stringNexOnOff,"N excess (ON - OFF) = %.d ", NexOnOff);
 sprintf(stringLiMaOnOff,"Signif (ON - OFF) = %.2f ", SignLiMaOnOff);


            bgoff2->SetParameters(aparOff, bparOff, normf/binwidth*binwidthOff);
            bgoff2->FixParameter(0, aparOff);
            bgoff2->FixParameter(1, bparOff);
            bgoff2->FixParameter(2, normf/binwidth*binwidthOff);
            bgoff2->SetLineColor(2);

   }

/*
      TCanvas canv("c1", "basic histograms", 600, 500);
      canv.SetBorderMode(0);
      canv.Divide(3,3);

     canv.cd(1);
     gPad->SetBorderMode(0);
     histlength.Draw();

     canv.cd(2);
     gPad->SetBorderMode(0);
     histwidth.Draw();
  
     canv.cd(3);
     gPad->SetBorderMode(0);
     gPad->SetLogx();
     gPad->SetLogy();
     histsize.Draw();

     canv.cd(4);
     gPad->SetBorderMode(0);
     histalpha.Draw();

     canv.cd(5);
     gPad->SetBorderMode(0);
     histdist.Draw();

     canv.cd(6);
     gPad->SetBorderMode(0);
     histmeanx.Draw();

     canv.cd(7);
     gPad->SetBorderMode(0);
     histmeany.Draw();

     canv.cd(8);
     gPad->SetBorderMode(0);
     hist2xy.Draw();

     canv.cd(9);
     gPad->SetBorderMode(0);
     histLoverS.Draw();

     canv.Modified();
     canv.Update();

     canv.DrawClone();




      TCanvas canvcor("c2", "correlation histograms", 600, 500);
      canvcor.SetBorderMode(0);
      canvcor.Divide(3,3);

     canvcor.cd(1);
     gPad->SetBorderMode(0);
     hist2lw.Draw();

     canvcor.cd(2);
     gPad->SetBorderMode(0);
     hist2lalpha.Draw();

     canvcor.cd(3);
     gPad->SetBorderMode(0);
     hist2ldist.Draw();

     canvcor.cd(4);
     gPad->SetBorderMode(0);
     hist2walpha.Draw();

     canvcor.cd(5);
     gPad->SetBorderMode(0);
     hist2wdist.Draw();

     canvcor.cd(6);
     gPad->SetBorderMode(0);
     hist2alphadist.Draw();

     canvcor.cd(7);
     gPad->SetBorderMode(0);
     gPad->SetLogy();
     hist2wsize.Draw();


     canvcor.cd(8);
     gPad->SetBorderMode(0);
     gPad->SetLogy();
     hist2alphasize.Draw();

     canvcor.cd(9);
     gPad->SetBorderMode(0);
     gPad->SetLogy();
     hist2distsize.Draw();


     canvcor.Modified();
     canvcor.Update();

     canvcor.DrawClone();
*/

/**********************************************************/
/* plot the alpha plot for the current sample */

    Char_t titelname[80];
    sprintf(titelname,"alpha plot. sample %d assumed source position: x = %.2f y = %.2f", num, XSOURCE, YSOURCE);


      TCanvas canval("c3", "canvas for alpha", 600, 500);
      canval.SetBorderMode(0);

     gPad->SetBorderMode(0);
     histalphafinal.SetMarkerStyle(20);
     histalphafinal.SetTitle(titelname);
     histalphafinal.SetFillColor(8);
     histalphafinal.Draw();
     bgoff->Draw("same");
     if (USEOFFDATA == kTRUE) bgoff2->Draw("same");

    leg = new TLegend(0.1,0.15,0.52,0.35);
//    leg->Draw();
    leg->AddEntry(fitbgpar,"fit for OFF region (30-90 deg)","l");
    leg->SetHeader("Legend");
    leg->SetFillColor(19);
//    leg->Draw();

    if (USEOFFDATA == kFALSE) 
    {
	    text = new TPaveText(0.53,0.45,0.9,0.65,"NDC");
	    text->AddText(0.4, 0.6, stringNex);
	    text->AddText(0.5, 0.3, stringlima);
	    text->SetTextSize(0.032);
	    text->Draw();
    }
    else
    {
            text = new TPaveText(0.53,0.25,0.9,0.55,"NDC");
	    text->AddText(0.4, 0.8, stringNex);
	    text->AddText(0.45, 0.6, stringlima);
            text->AddText(0.5, 0.4, stringNexOnOff);
            text->AddText(0.4, 0.2, stringLiMaOnOff);
            text->SetTextSize(0.032);
            text->Draw();
    }
     canval.Modified();
     canval.Update();

//     canval.DrawClone();

     TString strin = dirname + sourcename;
//     TString strin = "data/plots/alpha/dummySize";
     strin += SIZEMIN;
     strin += "Sample";
     if (num<10) strin +="0";
     strin += num;
     strin += ".root";


     canval.SaveAs(strin);  // please enable if you want to save alphaplots
cout << " alpha plot for the sample " << num << " has been saved into " << strin << endl;

     if(bgoff2) delete bgoff2;
     delete bgoff;

/**********************************************************/

/*
      TCanvas canvt("c4", "telescope", 600, 500);
      canvt.SetBorderMode(0);
      canvt.Divide(2,1);

     canvt.cd(1);
     gPad->SetBorderMode(0);
     histphi.Draw();

     canvt.cd(2);
     gPad->SetBorderMode(0);
     histtheta.Draw();

     canvt.DrawClone();

      TCanvas canvass("c5", "assymetry", 600, 500);
      canvass.SetBorderMode(0);
      canvass.Divide(2,1);

     canvass.cd(1);
     gPad->SetGridx();
     gPad->SetGridy();
     gPad->SetBorderMode(0);
     histAssym.Draw();

     canvass.cd(2);
     gPad->SetGridx();
     gPad->SetGridy();
     gPad->SetBorderMode(0);
     histAssymM3.Draw();

     canvass.DrawClone();

    TCanvas aftercutcanv("c1a", "basic histograms", 600, 500);
    aftercutcanv.SetBorderMode(0);
    aftercutcanv.Divide(3,3);

    aftercutcanv.cd(1);
    gPad->SetBorderMode(0);
    gPad->SetGridx();
    gPad->SetGridy();
    aftercuthistlength.Draw();

    aftercutcanv.cd(2);
    gPad->SetBorderMode(0);
    gPad->SetGridx();
    gPad->SetGridy();
    aftercuthistwidth.Draw();

    aftercutcanv.cd(3);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    gPad->SetLogy();
    gPad->SetGridx();
    gPad->SetGridy();
    aftercuthistsize.Draw();

    aftercutcanv.cd(4);
    gPad->SetBorderMode(0);
    gPad->SetGridx();
    gPad->SetGridy();
    aftercuthistalpha.Draw();

    aftercutcanv.cd(5);
    gPad->SetBorderMode(0);
    gPad->SetGridx();
    gPad->SetGridy();
    aftercuthistdist.Draw();

    aftercutcanv.cd(6);
    gPad->SetBorderMode(0);
    gPad->SetGridx();
    gPad->SetGridy();
    aftercuthistmeanx.Draw();

    aftercutcanv.cd(7);
    gPad->SetBorderMode(0);
    gPad->SetGridx();
    gPad->SetGridy();
    aftercuthistmeany.Draw();


    aftercutcanv.cd(8);
    gPad->SetBorderMode(0);
    gPad->SetGridx();
    gPad->SetGridy();
    aftercuthist2xy.Draw();

    aftercutcanv.Modified();
    aftercutcanv.Update();

    aftercutcanv.DrawClone();
*/
    Double_t rat1, rat2;
    rat1 = (double)AsGrNull / (double)AsLessNull;
    rat2 = (double)AsGrNullAfter / (double)AsLessNullAfter;

    cout << " Asymmetry M3 > 0 : " <<  AsGrNull << endl;
    cout << " Asymmetry M3 < 0 : " <<  AsLessNull << endl;
    cout << " Ratio (before cuts) : " <<  rat1 << endl;
    cout << " Asymmetry M3 > 0 (after): " <<  AsGrNullAfter << endl;
    cout << " Asymmetry M3 < 0 (after): " <<  AsLessNullAfter << endl;
    cout << " Ratio (after cuts) : " <<  rat2 << endl;


}

void callalphaplot()
{


   Int_t num;
   Double_t xpeakM, ypeakM, xpeakB, ypeakB;

   FILE *fp;
//   fp = fopen("data/trackMrk421_0505_2000.dat", "r");
   if (USEFILE == kTRUE) 
   {
	fp = fopen("data/trackCrab1502_1305_Berlin_5deg.dat", "r");

        for(Int_t i = 0; i < 2; i++)
        {
       		fscanf(fp,"%d %lf %lf %lf %lf", &num, &xpeakM, &ypeakM, &xpeakB, &ypeakB); 
       		cout << endl << " SUBS NUMBER " << num << ", xpeakM = " << xpeakM << 
                     ", ypeakM = " << ypeakM << endl;
       		cout << "                   xpeakB = " << xpeakB << ", ypeakB = " << ypeakB << endl;
       		if (num > 0 ) alphaplots(num, xpeakM, ypeakM);  
   	}

   	fclose(fp); 
cout << "FERTIG" << endl;
   }
   else  
     alphaplots(NUM, XSOUR, YSOUR);


/**********************************************************/
/* now calculate Nex and S for the overall alpha plot */:



    TF1 * fitbgparAll = new TF1("fbgparA", "[0]*x*x + [1]", 30., 90.);
    fitbgparAll->SetLineColor(2);
    fitbgparAll->SetLineWidth(3);

    histalphaAll.Fit("fbgparA","WR");

    Double_t apar = fitbgparAll->GetParameter(0);
    Double_t bpar = fitbgparAll->GetParameter(1);
    Double_t normf = 1.;

    TF1 * bgoff = new TF1("bgoff", parabfunc, 0., 90., 3);
    bgoff->SetParameters(apar, bpar, normf);
    bgoff->FixParameter(0, apar);
    bgoff->FixParameter(1, bpar);
    bgoff->FixParameter(2, normf);
    bgoff->SetLineColor(8);

// calc significance:
    Double_t Sign, Non, Noff, Nex;
    Double_t binwidth = histalphaAll.GetBinWidth(1);
    Double_t numbinMax = ALPHAMAX / binwidth;

    Non = 0.;
    for(Int_t i=1; i<=numbinMax;i++) Non += histalphaAll.GetBinContent(i);

//cout << histalphaAll.GetBinContent(1) + histalphaAll.GetBinContent(2) + 
//        histalphaAll.GetBinContent(3) + histalphaAll.GetBinContent(4) << endl;
    Noff = (1./3. * apar * pow(ALPHAMAX,3.) + bpar * ALPHAMAX)  /  binwidth;
    Nex = Non - Noff;

    Sign = LiMa17(Non,Noff,1.);

    cout << " Non : " << Non << " Noff : " << Noff << " Nex : " << Nex << endl;
    cout << " significance : " << Sign << " sigma" << endl;


    Char_t stringsig[80], stringNex[80];

    sprintf(stringsig,"Signif:  S = %.2f sigma", Sign);
    sprintf(stringNex,"N excess: Nex = %.d ", Nex);
/**********************************************************/


//  plot all alpha plots together

    TCanvas canvA("cA", "alphacanvas", 600, 500);
    canvA.SetBorderMode(0);

    gPad->SetBorderMode(0);
    histalphaAll.SetXTitle("alpha [deg]");
    histalphaAll.SetYTitle("Counts");
    histalphaAll.SetMarkerStyle(20);
    histalphaAll.SetFillColor(17);
    histalphaAll.Draw();
    bgoff->Draw("same");

//    leg = new TLegend(0.1,0.15,0.52,0.35);
//    leg->Draw();
//    leg->AddEntry(fitbgpar,"fit for OFF region (30-90 deg)","l");
//    leg->SetHeader("Legend");
//    leg->SetFillColor(19);
//    leg->Draw();

     text = new TPaveText(0.53,0.45,0.9,0.65,"NDC");
     text->AddText(0.4, 0.6, stringNex);
     text->AddText(0.45, 0.3, stringsig);
     text->SetTextSize(0.032);
     text->Draw();

     canvA.Modified();
     canvA.Update();


     Char_t xstr[20], ystr[20];
     if (USEFILE == kTRUE) 
     {
	sprintf(xstr,"%3.2f",0.);
	sprintf(ystr,"%3.2f",0.);
     }
     else
     {
	sprintf(xstr,"%3.2f",XSOUR);
	sprintf(ystr,"%3.2f",YSOUR);
     }


     TString string = dirname + sourcename;
     string += "X";
     string += xstr;
     string += "Y";
     string += ystr;
     string += "Size";
     string += SIZEMIN;
     string += "alpha";
     string += TMath::Nint(ALPHAMAX);
     string += ".root";

     canvA.SaveAs(string);
cout << " alpha plot has been saved into " << string << endl;

}