void ReadSetup(TString fname, MAstroCamera &cam)
{
    MMcConfigRunHeader *config=0;
    MGeomCam           *geom=0;

    TFile file(fname);
    TTree *tree = (TTree*)file.Get("RunHeaders");
    tree->SetBranchAddress("MMcConfigRunHeader", &config);
    if (tree->GetBranch("MGeomCam"))
        tree->SetBranchAddress("MGeomCam", &geom);
    tree->GetEntry(0);

    cam.SetMirrors(*config->GetMirrors());
    cam.SetGeom(*geom);
}

//--------------------------------------------------------------------------
//
//  findTelAxisFromStars
//
//  macro to 
//    - determine the positions of stars in the camera from the DC currents 
//    - and to find the expected position of the source in the camera
//

void findTelAxisFromStars(const TString filename="20040422_23213_D_Mrk421_E.root", const TString directory="/.magic/magicserv01/MAGIC/rootdata/2004_04_22/", const UInt_t numEvents = 0)
{
  gLog.SetNoColors();

  //-------------------------------------------------------------------
  // book some histograms which are to be filled in the event loop
  //

  TH2D *aberr = new TH2D("OpticalAberr", "Opt_Aberr",
                         100, 0.0, 650.0,  100, 0.0, 650.0);
  aberr->SetXTitle("Ideal distance from center [mm]");
  aberr->SetYTitle("Expected distance from center [mm]");

  TProfile *reldiff = new TProfile("RelativeDiff", "Rel_Diff",
				   100, 0.0, 650.0, -0.2, 0.2, "S");
  reldiff->SetXTitle("Ideal distance from center [mm]");
  reldiff->SetYTitle("(R - R0)/R0");

  //-------------------------------------------------------------------


  MParList  plist;
  MTaskList tlist;
  plist.AddToList(&tlist);

  MGeomCamMagic geomcam;
  MCameraDC     dccam;
  MStarCam      starcam;

  //$$$$$$$$$$$$$$$$ ww
  MObservatory  obs;
  MStarCam      sourcecam;
  sourcecam.SetName("MSourceCam");

  MHTelAxisFromStars htelaxis;

  plist.AddToList(&obs);
  plist.AddToList(&htelaxis);
  //$$$$$$$$$$$$$$$$ ww


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

  // Reads the trees of the root file and the analysed branches
  MReadReports read;
  read.AddTree("Currents"); 
  read.AddTree("Drive"); // If you do not include Drive info, the MFindStars class
                         // will not use the star catalog method
  read.AddFile(directory+filename); // after the reading of the trees!!!
  read.AddToBranchList("MReportCurrents.*");
  read.AddToBranchList("MReportDrive.*");

  MGeomApply geomapl;
  //  TString continuoslightfile =   
  //   "/data/MAGIC/Period016/rootdata/2004_04_16/20040416_22368_P_Off3c279-2CL100_E.root";

  TString continuoslightfile =
"/.magic/magicserv01/MAGIC/rootdata/2004_04_16/20040416_22368_P_Off3c279-2CL100_E.root";


  Float_t mindc = 0.9; //[uA]

  MCalibrateDC dccal;
  dccal.SetFileName(continuoslightfile);
  dccal.SetMinDCAllowed(mindc);

  const Int_t numblind = 5;
  const Short_t x[numblind] = { 47, 124, 470, 475, 571};
  const TArrayS blindpixels(numblind,(Short_t*)x);
  Float_t ringinterest = 100; //[mm]
  Float_t tailcut = 2.5;
  UInt_t integratedevents = 1;

  // We need the MAGIC mirror geometry from a MC header:
  //TString geometryfile = "/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root";
  TString geometryfile =
"/.magic/magicserv01/MAGIC/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root";




  // We need the Bright Star Catalog:
  //TString catalogfile = "/home/rwagner/bsc5.dat";
  TString catalogfile = "mtemp/mmpi/macros/bsc5.dat";

  //$$$$$$$$$$$$$$$$ ww
  MSourceDirections sdirs;
  sdirs.SetGeometryFile(geometryfile);

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 26);
  const Double_t dec = MAstro::Dms2Rad(38, 12, 36);
  sdirs.SetRaDec(ra,dec);
  sdirs.AddDirection(ra,dec,1,"Mkn 421");

  const Double_t ra  = MAstro::Hms2Rad(11,  0, 50);
  const Double_t dec = MAstro::Dms2Rad(39, 12, 44);
  sdirs.AddDirection(ra,dec,1,"My_UMa 49");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(38, 14, 29);
  sdirs.AddDirection(ra,dec,1,"My_UMa 51");

  // dummy source dirtections
  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(38, 20,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy1");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(38, 20,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy2");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(38, 30,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy3");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(38, 40,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy4");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(38, 50,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy5");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(39,  0,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy6");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(39, 10,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy7");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(39, 20,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy8");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(39, 30,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy9");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(39, 40,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy10");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(39, 50,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy11");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(40,  0,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy12");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(40, 10,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy13");

  const Double_t ra  = MAstro::Hms2Rad(11,  4, 31);
  const Double_t dec = MAstro::Dms2Rad(40, 20,  0);
  sdirs.AddDirection(ra,dec,1,"Dummy14");


  sdirs.SetRadiusFOV(3);
  //$$$$$$$$$$$$$$$$ ww

  MFindStars findstars;
  findstars.SetBlindPixels(blindpixels);
  findstars.SetRingInterest(ringinterest);
  findstars.SetDCTailCut(tailcut);
  findstars.SetNumIntegratedEvents(integratedevents);
  findstars.SetMinuitPrintOutLevel(-1);
  findstars.SetGeometryFile(geometryfile);
  findstars.SetBSCFile(catalogfile);
  const Double_t ra  = MAstro::Hms2Rad(11, 4, 26);
  const Double_t dec = MAstro::Dms2Rad(38, 12, 36);
  findstars.SetRaDec(ra,dec);
  findstars.SetLimMag(8);
  findstars.SetRadiusFOV(1.5);



  //$$$$$$$$$$$$$$$$ ww
  Int_t InputType = 0;
  //findstars.SetUseCorrelatedGauss(kFALSE);

  MTelAxisFromStars telaxis;
  //telaxis.FixRotationAngleAt(-1.0);
  //telaxis.FixScaleFactorAt(-1.0);
  telaxis.SetInputType(InputType);

  MFillH fillhisto("MHTelAxisFromStars[MHTelAxisFromStars]","");
  htelaxis.SetInputType(InputType);
  //$$$$$$$$$$$$$$$$ ww


  tlist.AddToList(&geomapl);
  tlist.AddToList(&read);
  tlist.AddToList(&dccal);
  tlist.AddToList(&findstars, "Currents");

  //$$$$$$$$$$$$$$$$ ww
  tlist.AddToList(&sdirs);
  tlist.AddToList(&telaxis);
  tlist.AddToList(&fillhisto);
  //$$$$$$$$$$$$$$$$ ww


  
  // The following lines you only need if in addition you want to display
  // independent MAstroCamera output
  //
  //  TString fname = "/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root";
  //  MObservatory magic1;      
  //  const Double_t ra  = MAstro::Hms2Rad(11, 4, 26); //Mkn421
  //  const Double_t dec = MAstro::Dms2Rad(38, 12, 36);
  //
  //  MAstroCamera stars;        
  //  ReadSetup(fname, stars); 
  //  stars.SetLimMag(9);
  //  stars.SetRadiusFOV(3);
  //  stars.SetRaDec(ra, dec);
  //  stars.ReadBSC("/home/rwagner/bsc5.dat");
  //  stars.SetObservatory(magic1);

  MEvtLoop evtloop;
  evtloop.SetParList(&plist);
     
  if (!evtloop.PreProcess())
    return;
  
  MHCamera display(geomcam);
  display.SetPrettyPalette();
  display.Draw();
  gPad->cd(1);
  starcam.Draw();
  
  UInt_t numevents=0;
  
  while (tlist.Process())
    {
      //gLog << "---------------------------------------------" << endl;
      //gLog << "Macro : MStarPos content = " << endl;
      //starcam.Print("namepossizchierr");
      //gLog << "---------------------------------------------" << endl;
      //gLog << "Macro : MSourcePos content = " << endl;
      //sourcecam.Print("namepossizchierr");
      //gLog << "---------------------------------------------" << endl;

      numevents++;
      if (numevents%integratedevents==0)
	{
	  display.SetCamContent(findstars.GetDisplay());
	  gPad->Modified();
	  gPad->Update();	      
        // This line prints the results:
	// 	  starcam.Print();

	// This is how to access the TList of stars:
	//	  TList* starlist = starcam.GetList();

	// This is how to iterate over stars found:
	//	  TIter Next(starlist);
	// 	  MStarPos* star;
	// 	  UInt_t starnum = 0;
	//        cout << filename << " ";
	//        cout << "Iterating over list" << endl;
	// 	  while ( (star=(MStarPos*)Next()) ) 
	// 	  {
	//          cout << "star[" << starnum << "] ";
	//          cout << star->GetMeanX() << " " 
	//               << star->GetMeanY() << " ";
	//	    starnum++;
	//	  }
	//        cout << endl;
	}

      // fill the histograms
      // loop over stars in "MSourceCam"
      TList* sourcelist = sourcecam.GetList();
      TIter Nextsource(sourcelist);
      MStarPos* source;
      while ( (source=(MStarPos*)Nextsource()) ) 
      {
        Double_t R = sqrt( source->GetXExp()*source->GetXExp()
      		   + source->GetYExp()*source->GetYExp() );

        Double_t R0= sqrt( source->GetXIdeal()*source->GetXIdeal()
      		   + source->GetYIdeal()*source->GetYIdeal() );

        //cout << "R0, R = " << R0 << ",  " << R << endl;

	aberr->Fill(R0, R, 1.0);
        reldiff->Fill(R0, (R-R0)/R0), 1.0);
      }
	
      //--------------------------------------
      MTime time;
      time.Set(2004, 4, 22, 21, 51, 15);
      
      //superimpose star picture
      //       stars.SetTime(time);
      //       TObject *o = stars.Clone();
      //       o->SetBit(kCanDelete);
      //       o->Draw();
      
      // wait after each event
      //if (!HandleInput())
      //  break;

    }



  evtloop.PostProcess();
  tlist.PrintStatistics();
    
  plist.Print();

  //$$$$$$$$$$$$$$$$ ww
  
  //gLog << "Event loop finished; call DrawClone of MHTelAxisFromStars" << endl;

  //TObject *srccam = plist.FindObject("MSourceCam");
  //gLog << "srccam = " << srccam << endl;

  TObject *obj = plist.FindObject("MHTelAxisFromStars");
  if (obj)
  {
    //obj->Print();
    //obj->Dump();
    obj->DrawClone();
  }
  else
    gLog << "address of MHTelAxisFromStars container is zero" << endl;
  
  //$$$$$$$$$$$$$$$$ ww

  // plot the histograms
  TCanvas *canv = new TCanvas("Optical_Aberration", "Optical aberration",
			      0, 0, 600, 300);
  canv->Divide(2,1);
  canv->SetBorderMode(0); 
  
  canv->cd(1);
  aberr->Draw();

  canv->cd(2);
  reldiff->Draw();

}

//----------------------------------------------------------------------
//
//

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;
}
//=========================================================================









