source: trunk/MagicSoft/Mars/mtemp/MFindStars.cc@ 3908

/* ======================================================================== *\
!
! *
! * 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): Robert Wagner, 2/2004 <mailto:rwagner@mppmu.mpg.de>
!   Author(s): Javier López , 4/2004 <mailto:jlopez@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MFindStars
//
/////////////////////////////////////////////////////////////////////////////
#include "MFindStars.h"

#include <TFile.h>
#include <TTree.h>

#include "MObservatory.h"
#include "MAstroCamera.h"
#include "MMcConfigRunHeader.h"

#include "MLog.h"
#include "MLogManip.h"

#include "MHCamera.h"

#include "MGeomCam.h"
#include "MGeomPix.h"
#include "MCameraDC.h"
#include "MTime.h"
#include "MReportDrive.h"
#include "MStarLocalCam.h"
#include "MStarLocalPos.h"

#include "MParList.h"

ClassImp(MFindStars);
using namespace std;

MFindStars::MFindStars(const char *name, const char *title)
{
  fName  = name  ? name  : "MFindStars";
  fTitle = title ? title : "Tool to find stars from Base Currents";

  fMaxNumIntegratedEvents = 10;
  fRingInterest = 100; //[mm] ~ 0.3 deg

  fPixelsUsed.Set(577);
  fPixelsUsed.Reset((Char_t)kTRUE);
}

Int_t MFindStars::PreProcess(MParList *pList)
{

    fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));

    if (!fGeomCam)
    {
      *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
      return kFALSE;
    }

    fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));

    if (!fCurr)
    {
      *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
      return kFALSE;
    }

    fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents"));

    if (!fTimeCurr)
    {
      *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl;
      return kFALSE;
    }

   fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));

    if (!fDrive)
    {
      *fLog << err << AddSerialNumber("MReportDrive") << " not found ... aborting" << endl;
      return kFALSE;
    }

    fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam"));
    if (!fStars)
    {
      *fLog << err << AddSerialNumber("MStarLocalCam") << " cannot be created ... aborting" << endl;
      return kFALSE;
    }

//Initialitation MAstroCamera
    // Name of a MC file having MGeomCam and MMcConfigRunHeader
    TString fname = "../data/Gamma_zbin9_90_7_1480to1489_w0.root";

    // Time for which to get the picture
//    MTime time;
//    time.Set(2004, 2, 28, 01, 32, 15);

    // Current observatory
    MObservatory magic1;

    // Right Ascension [h] and declination [deg] of source
    // Currently 'perfect' pointing is assumed
//    const Double_t ra  = MAstro::Hms2Rad(5, 34, 31.9);
//    const Double_t dec = MAstro::Dms2Rad(22, 0, 52.0);

      // --------------------------------------------------------------------------
    // Create camera display from geometry
    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);

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


    fAstro.SetLimMag(6);
    fAstro.SetRadiusFOV(3);
    fAstro.ReadBSC("../data/bsc5.dat");

    fAstro.SetObservatory(magic1);
    // Time for which to get the picture
//    MTime time;
//    time.Set(2004, 2, 28, 01, 32, 15);
//   fAstro.SetTime(time);
    fAstro.SetTime(*fTimeCurr);
    fAstro.SetGuiActive();

    fAstro.SetStarList(fStars->GetList());

    return kTRUE;
}

Int_t MFindStars::Process()
{

    Float_t ra  = fDrive->GetRa();
    Float_t dec = fDrive->GetDec();

    fAstro.SetRaDec(ra, dec);
    fAstro.SetGuiActive();
   
    fAstro.FillStarList();

    if (fNumIntegratedEvents < fMaxNumIntegratedEvents)
    {
      fDisplay.AddCamContent(*fCurr);
      fNumIntegratedEvents++;
    }
    else
    {
        //loop to extract position of stars on the camera
        TIter Next(fStars->GetList());
        MStarLocalPos* star;
        while ((star=(MStarLocalPos*)Next())) 
            FindStar(star);

        //After finding stars reset all vairables
        fDisplay.Reset();
        fNumIntegratedEvents=0;
    }

  return kTRUE;
}

Int_t MFindStars::PostProcess()
{
    fStars->Print();
    return kTRUE;
}

void MFindStars::SetBlindPixels(TArrayS blindpixels)
{
    Int_t npix = blindpixels.GetSize();

    for (Int_t idx=0; idx<npix; idx++)
        fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;

    fDisplay.SetUsed(fPixelsUsed);
}

Bool_t MFindStars::FindStar(MStarLocalPos* star)
{    

    Float_t expX = star->GetXExp();
    Float_t expY = star->GetYExp();

    UInt_t numPixels = fGeomCam->GetNumPixels();

// First define a area of interest around the expected position of the star
    for (UInt_t pix=1; pix<numPixels; pix++)
    {

        Float_t pixXpos=(*fGeomCam)[pix].GetX();
        Float_t pixYpos=(*fGeomCam)[pix].GetY();
        Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
                            (pixYpos-expY)*(pixYpos-expY));

      if (dist < fRingInterest)
          fPixelsUsed[pix]=(Char_t)kTRUE;
    }

    fDisplay.SetUsed(fPixelsUsed);

// Find the two close pixels with the maximun dc
    for (UInt_t pix=1; pix<numPixels; pix++)
    {
        Float_t maxDC = 0;

        if(fDisplay.IsUsed(pix))
        {
            Float_t dc[2];
            Float_t maxPix[2];
            dc[0] = fDisplay.GetBinContent(pix+1);
            
            MGeomPix g = (*fGeomCam)[pix];
            Int_t numNextNeighbors = g.GetNumNeighbors();
            
            Float_t dcsum;
            for(Int_t nextNeighbor=0; nextNeighbor<numNextNeighbors; nextNeighbor++)
            {
                if(fDisplay.IsUsed(pix))
                {
                    UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
                    dc[1] = fDisplay.GetBinContent(swneighbor+1);
                    
                    
                    dcsum = dc[0] + dc[1];
                    
                    if(dcsum > maxDC*2)
                    {
                        maxPix[0] = pix;
                        maxPix[1] = swneighbor;
                        maxDC = dcsum/2;
                    }   
                }
            }
        }
    }

    // determine mean x and mean y of the selected px
    Float_t meanX=0;
    Float_t meanY=0;
    Float_t meanSqX=0;
    Float_t meanSqY=0;
    Float_t sumCharge=0;
    UInt_t usedPx=0;    
    for(UInt_t pix=0; pix<numPixels; pix++)
    {
        if(fDisplay.IsUsed(pix))
        {
            usedPx++;

            Float_t charge=fDisplay.GetBinContent(pix+1);
            Float_t pixXpos=(*fGeomCam)[pix].GetX();
            Float_t pixYpos=(*fGeomCam)[pix].GetY();

            meanX+=charge*pixXpos;
            meanY+=charge*pixYpos;
            meanSqX+=charge*pixXpos*pixXpos;
            meanSqY+=charge*pixYpos*pixYpos;
            sumCharge+=charge;
//          fDisplay.ResetUsed(i); // this px must not be used again!
        } //for... use  
    }

    star->SetCalcValues(sumCharge,meanX,meanY,meanSqX,meanSqY);
    fStars->GetList()->Add(star);
    
    return kTRUE;
}