source: trunk/MagicSoft/Mars/mtemp/meth/MDCA.cc@ 4433

/* ======================================================================== *\
!
! *
! * 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): Sebastian Commichau 05/2004 <mailto:commichau@particle.phys.ethz.ch>
!   Author(s): Sabrina Stark       05/2004 <mailto:lstark@particle.phys.ethz.ch>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

// Container to store the DCA stuff - it offers a nice draw option...                   


#include "MDCA.h"

using namespace std;

ClassImp(MDCA); 

// Default constructor
MDCA::MDCA(const char *name, const char *title) : fXRef(0.0), fYRef(0.0)
{
    fName = name ? name : "MDCA";
    fTitle = title ? title : "Storage container for Hillas parameters and DCA of one event";
 
    Reset();
   
    fEllipse = new TEllipse;
    fRefCircle = new TEllipse;
    
    fLineL = new TLine;
    fLineW = new TLine;
    fLineX = new TLine;
    fLineY = new TLine;
    fLineDCA = new TLine;
    fLineMean = new TLine;
  
}

// Destructor: Deletes ellipse and lines if they do exist

MDCA::~MDCA()
{
    Clear();
}

// Initialize parameters with default values

void MDCA::Reset()
{

    fLength = -1;
    fWidth  = -1;
    fDelta0 =  0;
    fMeanX  =  0;
    fMeanY  =  0;
    fDelta1 =  0;
    fDCA    = -1;
    fX1W    =  0;
    fY1W    =  0;
    fX2W    =  0;
    fY2W    =  0;
    fX1L    =  0;
    fY1L    =  0;
    fX2L    =  0;
    fY2L    =  0;
    fXDCA   =  0;
    fYDCA   =  0;
  
}


// Print parameters to *fLog

void MDCA::Print(Option_t *) const
{
    Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg;

    if (atg<0)
        atg += 180;

    *fLog << all;
    *fLog << "Basic Image Parameters (" << GetName() << ")" << endl;
    *fLog << " - Length      [mm]  = " << fLength << endl;
    *fLog << " - Width       [mm]  = " << fWidth  << endl;
    *fLog << " - Delta0      [deg] = " << fDelta0*kRad2Deg << endl;
    *fLog << " - Meanx       [mm]  = " << fMeanX  << endl;
    *fLog << " - Meany       [mm]  = " << fMeanY  << endl;
    *fLog << " - atg(y/x)    [deg] = " << atg     << endl;
    *fLog << " - DCA         [mm]  = " << fDCA << endl;
    *fLog << " - Delta1      [deg] = " << fDelta1*kRad2Deg << endl;

    
}

void MDCA::Paint(Option_t *opt)
{
    Clear();

    if (fLength<=0 || fWidth<=0) //fLength<0 || fWidth<0 doesn't look nice...
        return;                  //We get a circle!

    // Length line
    fLineL = new TLine(fX1L, fY1L, fX2L, fY2L);
    fLineL->SetLineWidth(2);
    fLineL->SetLineColor(2);
    fLineL->Draw();

    // Width line
    fLineW = new TLine(fX1W, fY1W, fX2W, fY2W);
    fLineW->SetLineWidth(2);
    fLineW->SetLineColor(2);
    fLineW->Draw();

    // Coordinate system
    fLineX = new TLine(-600,fYRef,600,fYRef);
    fLineY = new TLine(fXRef,-600,fXRef,600);
    fLineX->SetLineWidth(2);
    fLineX->SetLineColor(1);
    fLineY->SetLineWidth(2);
    fLineY->SetLineColor(1);
    fLineX->Draw();
    fLineY->Draw();
    
    // DCA line
    fLineDCA = new TLine(fXRef,fYRef,fXDCA+fXRef,fYDCA+fYRef);
    fLineDCA->SetLineWidth(2);
    fLineDCA->SetLineColor(2);
    fLineDCA->Draw();

    // COG line
    fLineMean = new TLine(fXRef,fYRef,fMeanX,fMeanY);
    fLineMean->SetLineWidth(2);
    fLineMean->SetLineColor(2);
    fLineMean->Draw();

    // Reference point marker
    fRefCircle = new TEllipse(fXRef, fYRef, 5, 5, 0, 360, 0);
    fRefCircle->SetLineColor(8);
    fRefCircle->SetFillColor(8);
    fRefCircle->Draw();

    // Hillas ellipse
    fEllipse = new TEllipse(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta0*kRad2Deg+180);
    fEllipse->SetLineWidth(2);
    fEllipse->SetLineColor(2);
    fEllipse->Draw();

}


// If an ellipse and lines exist they will be deleted 

void MDCA::Clear(Option_t *)
{
    if (!fEllipse && !fRefCircle && !fLineL && !fLineW && !fLineX && !fLineY && !fLineDCA && !fLineMean) 
        return;
 
    delete fEllipse;
    delete fRefCircle;
    delete fLineL;
    delete fLineW;
    delete fLineX;
    delete fLineY;
    delete fLineDCA;
    delete fLineMean;
     
    fLineL = NULL;
    fLineX = NULL;
    fLineY = NULL;
    fLineW = NULL;
    fLineDCA = NULL;
    fLineMean = NULL;
    fEllipse = NULL;
    fRefCircle = NULL;
}


Int_t MDCA::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hil)
{
    // Get basic Hillas parameters from MHillas
    fDelta0 = hil.GetDelta();
    fMeanX  = hil.GetMeanX();
    fMeanY  = hil.GetMeanY();
    fLength = hil.GetLength();
    fWidth  = hil.GetWidth();

    // The Length Line - rotation and shift
    fX1L = - (fLength+OffsetL)*cos(fDelta0) + fMeanX; // [mm]    
    fY1L = - (fLength+OffsetL)*sin(fDelta0) + fMeanY; // [mm]    
    fX2L =   (fLength+OffsetL)*cos(fDelta0) + fMeanX; // [mm]    
    fY2L =   (fLength+OffsetL)*sin(fDelta0) + fMeanY; // [mm]

    // The Width Line - rotation and shift
    fX1W =   (fWidth+OffsetW)*sin(fDelta0) + fMeanX;  // [mm]   
    fY1W = - (fWidth+OffsetW)*cos(fDelta0) + fMeanY;  // [mm]   
    fX2W = - (fWidth+OffsetW)*sin(fDelta0) + fMeanX;  // [mm]   
    fY2W =   (fWidth+OffsetW)*cos(fDelta0) + fMeanY;  // [mm]  

    // Vector of orientation of the shower axis
    fr1 = fX2L-fX1L;  
    fr2 = fY2L-fY1L;  
       
    // Determine parameters to calculate coordinates of the DCA vector
    flambda = (fr1*(fXRef-fMeanX) + fr2*(fYRef-fMeanY))/(fr1*fr1 + fr2*fr2);
    fmu = (fMeanY-fYRef)/fr1 + flambda*fr2/fr1;

    // Components of the DCA vector
    fXDCA = -fmu*fr2;
    fYDCA =  fmu*fr1;
  
    // Components of vector going from intersection point of the DCA vector
    // with the shower axis to the COG
    fd1 = fMeanX + fmu*fr2 - fXRef;
    fd2 = fMeanY - fmu*fr1 - fYRef;

    // Calculate DCA value 
    fDCA = sqrt(fXDCA*fXDCA + fYDCA*fYDCA);

    // Calculate angle of the shower axis with respect to the x-axis
    fDelta1 = acos(fd1/sqrt(fd1*fd1 + fd2*fd2));

    // Calculate angle of the shower axis with respect to the y-axis
    //fDelta1 = acos(fd2/sqrt(fd1*fd1 + fd2*fd2)); 

    // Determine the correct sign of the DCA (cross product of DCA vector and 
    // vector going from the intersection point of the DCA vector with the shower axis
    // to the COG)
    if((fmu*(-fr2*(fMeanY-fYRef)-fr1*(fMeanX-fXRef)))<0)
       fDCA = -fDCA;    

    gRandom->Rannor(gx,gy);
    gx = fabs(gx);
  
    // This is nice but does not remove the systematics in the profile plot...
    //if(((1-0.6*gx)*(180-kRad2Deg*fDelta1)>120) || ((1-0.6*gx)*kRad2Deg*fDelta1>120))
      // fDCA = -1;

    // Enlarge the interval of Delta correctly...
    if((fMeanY-fYRef-fmu*fr1)<0)
       fDelta1 = TwoPi-fDelta1;

    // Enlarge the interval of Delta correctly... (Delta with respect to the y-axis)
    // if(-(fMeanX-fXRef+fmu*fr2)<0)
     //  fDelta1 = TwoPi-fDelta1;

    // This has to be improved...
    if(fr1 == 0 || fr2 == 0 || (fr1*fr1+fr2*fr2) == 0 || sqrt(fd1*fd1+fd2*fd2) == 0)
       fDCA = -fDCA;

    SetReadyToSave();
   
    return 0;
}

void MDCA::SetRefPoint(const Float_t fXRef0, const Float_t fYRef0)
{
    fXRef = fXRef0;
    fYRef = fYRef0;
}