/* ======================================================================== *\
!
! *
! * 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 Cerenkov 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): Luisa Sabrina Stark Schneebeli, 07/2004 <stark@particle.phys.ethz.ch>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!   This macro is partially derived from findstars.C
\* ======================================================================== */


Bool_t HandleInput()
{
    TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
    while (1)
    {
        //
        // While reading the input process gui events asynchronously
        //
        timer.TurnOn();
        TString input = Getline("Type 'q' to exit, <return> to go on: ");
        timer.TurnOff();

        if (input=="q\n")
            return kFALSE;

        if (input=="\n")
            return kTRUE;
    };

    return kFALSE;
}


void CorrectSrcPos(const TString filename="dc_2004_04_21_22_*_Mrk421.root", const TString directory="/scratch/Period016/cacodata/2004_04_22/", const UInt_t numEvents = 0)
{

  //
  // Create a empty Parameter List and an empty Task List
  // The tasklist is identified in the eventloop by its name
  //
  MParList  plist;
  
  MTaskList tlist;
  plist.AddToList(&tlist);


  MGeomCamMagic geomcam;
  MCameraDC     dccam;
  MStarLocalCam starcam;

  plist.AddToList(&geomcam);
  plist.AddToList(&dccam);
  plist.AddToList(&starcam);

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

  // Reads the trees of the root file and the analysed branches
  MReadReports read;
  read.AddTree("Currents"); 
  read.AddFile(directory+filename);     // after the reading of the trees!!!
  read.AddToBranchList("MReportCurrents.*");

  MGeomApply geomapl;
  TString continuoslightfile = 
    "/scratch/Period016/cacodata/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
  MCalibrateDC dccal;
  dccal.SetFileName(continuoslightfile);
  

  const Int_t numblind = 12;
  const Short_t x[numblind] = {  8,  27,
                               507, 508, 509, 510, 511,
                               543, 559, 560, 561, 567};
  const TArrayS blindpixels(numblind,(Short_t*)x);
  Float_t ringinterest = 100; //[mm]
  Float_t tailcut = 3.5;
  UInt_t integratedevents = 5;

  MFindStars findstars;
  //  findstars.SetBlindPixels(blindpixels);
  findstars.SetRingInterest(ringinterest);
  findstars.SetDCTailCut(tailcut);
  findstars.SetNumIntegratedEvents(integratedevents);
  findstars.SetMinuitPrintOutLevel(0);

  findstars.SetSourceRaH(11);   //  Mrk421
  findstars.SetSourceRaM(4);
  findstars.SetSourceRaS(26.0);
  findstars.SetSourceDecD(38);
  findstars.SetSourceDecM(12);
  findstars.SetSourceDecS(36.0);


  // prints
  MPrint pdc("MCameraDC");
  MPrint pstar("MStarLocalCam");
  
  tlist.AddToList(&geomapl);
  tlist.AddToList(&read);
  tlist.AddToList(&dccal);
  //  tlist.AddToList(&pdc, "Currents");
  tlist.AddToList(&findstars, "Currents");
  //  tlist.AddToList(&pstar, "Currents");

  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop;
  evtloop.SetParList(&plist);
     
  const Int_t numpar = 3;
  Float_t starpos[10*numpar];       
  Float_t starposcat[10*numpar];       
  
  if (numEvents > 0)
  {
      if (!evtloop.Eventloop(numEvents))
	  return;
  }
  else
  {
      if (!evtloop.PreProcess())
	  return;
      
       MHCamera display(geomcam);
       display.SetPrettyPalette();
       
       TCanvas *c = new TCanvas("c", "FindStars", 100, 100,600,600); 
       
       TH2D *h0 = new TH2D("h0","",1,-2,2,1,-2,2);
       h0->GetXaxis()->SetTitle("(deg)");
       h0->GetYaxis()->SetTitle("(deg)");
       h0->Draw();
       
       display.Draw();
       gPad->cd(1);
       //display.Draw("same");
       starcam.Draw();

       UShort_t year;
       Byte_t month, day;

       MTime *evttime = (MTime*)plist->FindObject("MTimeCurrents");

       TList *starcatalog = new TList;
       TList *stardc = new TList;
       TVector3 *res = new TVector3;

       Int_t maxdist = 100;
       UInt_t numevents=0;
       Int_t numEvt=0;
       const Int_t maxnumEvt = 10000;
       Int_t col = 2;
       const Double_t kDeg2mm = 189/0.6;

       Double_t resx[maxnumEvt];
       Double_t resy[maxnumEvt];
       Double_t time[maxnumEvt];

       
//
// Find in the file name the date, run and project name 

       Size_t pos = filename.Last('/');
       TString iRun = TString(filename(pos+24,5));
       TString iYear = TString(filename(pos+4,4));
       TString iMonth = TString(filename(pos+9,2));
       TString iDay = TString(filename(pos+12,2));
       
       TString iHour = TString(filename(pos+15,2));
       TString iMin = TString(filename(pos+18,2));
       TString iSec = TString(filename(pos+21,2));
       
       Size_t poslast = filename.Last('.');
       Size_t posfirst = poslast-1;
       while (filename[posfirst] != '_')
	 posfirst--;
       
       TString iSource = TString(filename(posfirst+1,poslast-posfirst-1));
       char str[100];
       char str2[100];
       sprintf(str,"Date %s %s %s  Run %s   Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data());
       sprintf(str2,"CorrectSrcPos_%s%s%s_%s_%s.txt",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data());
       
       UInt_t evtyear, evtmonth, evtday, evthour, evtmin, evtsec;
       UShort_t evtmsec;

       ofstream fposout1;
       //       fposout1.open (str2, ofstream::app);
       fposout1.open (str2, ofstream::trunc);

       while (tlist.Process())
         {
           if (numevents >= integratedevents)
             {
	       evttime->GetTime(evthour, evtmin, evtsec,evtmsec);
	       evttime->Print();
	       evthour = evttime->Hour();
	       evtmin = evttime->Min();
	       evtsec = evttime->Sec();

               //               display.SetCamContent(dccam);
               display.SetCamContent(findstars.GetDisplay());
	       gPad->Modified();
	       gPad->Update();
	       Int_t numStars = starcam.GetNumStars();
	       starcam.Print("bkmaxposrotsizechi");

	       res = findstars.GetCorrSrcPos();
	       cout<<"Corrected source position = "<<res->X()<<" , "<<res->Y()<<endl;
	       resx[numEvt] = res->X();
	       resy[numEvt] = res->Y();
	       time[numEvt] = numEvt;
	       numEvt ++;

	       // write the corrected source position and the star positions to file
	       if (!fposout1)
		 cout << "Could not open output file to write the center of the fit." << endl;
	       else{
		 for ( Int_t i = 0; i < numStars; i++)
		   {
		     if (i == 0){
		       fposout1 << str<<" "<<evthour<< " "<<evtmin<<" "<<evtsec<<" "<<evtmsec;
		       fposout1 <<" "<< res->X()<< " " << res->Y();
		     }
		     fposout1 <<" "<< starcam[i].GetMeanX()<< " " << starcam[i].GetMeanY();
		   }
		 }
	       fposout1 << endl;
	       numevents = 0;

               // Remove the comments if you want to go through the file
               // event-by-event:
	       //if (!HandleInput())
	       // break;
             }
           numevents++;
	 } 
       fposout1.close();

       // Print source position
       gROOT->Reset();
       gStyle->SetCanvasColor(0);
       gStyle->SetCanvasBorderMode(0);
       gStyle->SetPadBorderMode(0);
       gStyle->SetFrameBorderMode(0);
       gStyle->SetOptStat(00000000);

       TCanvas *c1 = new TCanvas("c1", "SourcePosition", 200, 200,600,1200); 
       c1->Divide(1,3);

       TMath math;
    
       Double_t minresx, maxresx;
       minresx = resx[math.LocMin(numEvt,resx)];
       maxresx = resx[math.LocMax(numEvt,resx)];
       
       Double_t minresy, maxresy;
       minresy = resy[math.LocMin(numEvt,resy)];
       maxresy = resy[math.LocMax(numEvt,resy)];
       
       Double_t minres, maxres;
       minres = math.Min(minresx,minresy);
       maxres = math.Max(maxresx,maxresy);

       Double_t diff;
       diff = maxres - minres;
       diff = 0.1*diff;
       minres = minres - diff;
       maxres = maxres + diff;

       Double_t mintime, maxtime;
       mintime = time[math.LocMin(numEvt,time)];
       maxtime = time[math.LocMax(numEvt,time)];

       c1->cd(1);
       TH2D *h1 = new TH2D("h1","",1,mintime-1,maxtime+1,1,minres,maxres);
       h1->GetXaxis()->SetTitle("(evt.nr)");
       h1->GetYaxis()->SetTitle("(mm)");
       h1->Draw();

       TGraph *grtimeevolx = new TGraph(numEvt,time,resx);
       grtimeevolx->SetMarkerColor(215);
       grtimeevolx->SetMarkerStyle(20);
       grtimeevolx->SetMarkerSize (0.4);
       grtimeevolx->Draw("P");
       
       TGraph *grtimeevoly = new TGraph(numEvt,time,resy);
       grtimeevoly->SetMarkerColor(6);
       grtimeevoly->SetMarkerStyle(24);
       grtimeevoly->SetMarkerSize (0.4);
       grtimeevoly->Draw("P");
       
       legxy = new TLegend(0.8,0.85,0.95,0.95);
       legxy.SetTextSize(0.03);
       legxy.AddEntry(grtimeevolx,"src pos x","P");
       legxy.AddEntry(grtimeevoly,"src pos y","P");
       legxy.Draw();

       c1->cd(2);
       TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minresx,maxresx);
       h2->GetXaxis()->SetTitle("(evt.nr)");
       h2->GetYaxis()->SetTitle("(mm)");
       h2->Draw();

       TGraph *grtimeevolx = new TGraph(numEvt,time,resx);
       grtimeevolx->SetMarkerColor(215);
       grtimeevolx->SetMarkerStyle(20);
       grtimeevolx->SetMarkerSize (0.4);
       grtimeevolx->Draw("P");
       
       legxy = new TLegend(0.8,0.85,0.95,0.95);
       legxy.SetTextSize(0.03);
       legxy.AddEntry(grtimeevolx,"src pos x","P");
       legxy.Draw();

       c1->cd(3);
       TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minresy,maxresy);
       h3->GetXaxis()->SetTitle("(evt.nr)");
       h3->GetYaxis()->SetTitle("(mm)");
       h3->Draw();

       TGraph *grtimeevoly = new TGraph(numEvt,time,resy);
       grtimeevoly->SetMarkerColor(6);
       grtimeevoly->SetMarkerStyle(24);
       grtimeevoly->SetMarkerSize (0.4);
       grtimeevoly->Draw("P");
       
       legxy = new TLegend(0.8,0.85,0.95,0.95);
       legxy.SetTextSize(0.03);
       legxy.AddEntry(grtimeevoly,"src pos y","P");
       legxy.Draw();

       evtloop.PostProcess();
  }

  tlist.PrintStatistics();

}


