source: trunk/MagicSoft/Mars/mtemp/mifae/programs/srcPos.cc@ 4019

////////////////////////////////////////////////////////////////////////////////////
//
//             _____ Source Position macro_____
//
//  Take as input root files with hillas parameters and recompute the ones depending
//  on the source position for a new (input) position, optionally rotating it in an
//  event by event basis. Output is a file with recomputed hillas parameters
//
//                 Ester Aliu   <aliu@ifae.es>
//                 Oscar Blanch <blanch@ifae.es>
//                 Javier Rico  <jrico@ifae.es>
////////////////////////////////////////////////////////////////////////////////////

#include <fstream>
#include <iostream>

#include "TString.h"

#include "MHillasSrcCalc.h"
#include "MHillasSrc.h"
#include "MSrcRotate.h"
#include "MSrcTranslate.h"
#include "MParList.h"
#include "MTaskList.h"
#include "MHillas.h"
#include "MReadTree.h"
#include "MEvtLoop.h"
#include "MLog.h"
#include "MArgs.h"
#include "MWriteRootFile.h"

using namespace std;

Bool_t readDatacards(TString& filename);
void srcPos();

//-----------------------------------------------------------------------------
// declaration of variables read from datacards
//-----------------------------------------------------------------------------

TString  inputFile;
TString  outputFile;
ULong_t  nmaxevents=999999999;
Float_t  xsrcpos=0.;
Float_t  ysrcpos=0.;
Bool_t   kRotate=1;
Bool_t   kSrcPolicy=kTRUE;
Double_t fRA= -1.;
Double_t fDEC= -1.;

//-----------------------------------------------------------------------------
// constants
//-----------------------------------------------------------------------------

const TString defaultcard="srcpos.datacard";

const Float_t alphamin = 0.;       // minimum value in alpha (degrees) 
const Float_t alphamax = 90.;      // maximum value in alpha (degrees)
const Float_t conver   = 189./0.6; // conversion factor degrees to mm

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

static void Usage()
{
  gLog <<endl;
  gLog << "Usage is:" << endl;
  gLog << "   srcPos [-h] [-?] <datacards>" << endl << endl;
  gLog << "     <datacards>: datacards file name (dafault " << defaultcard <<")" << endl;
  gLog << "     -?/-h: This help" << endl << endl;
}

//-----------------------------------------------------------------------------
int main(int argc, char **argv)
{
  // evaluate arguments
  MArgs arg(argc, argv);
  if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
    {
      Usage();
      return -1;
    }

  TString datacard  = arg.GetArgumentStr(0);
  if(!datacard.Length())
    datacard = defaultcard;

  if(!readDatacards(datacard))
    {
      cout << "Error reading datacards. Stoping" << endl;
      return -1;
    }
  srcPos();
}

//-----------------------------------------------------------------------------
void srcPos()
{
  // variable declaration  
  Float_t xpos=xsrcpos*conver; // [mm]
  Float_t ypos=ysrcpos*conver; // [mm] 
  
  // containers
  MParList  plist;  
  MTaskList tlist;

  // include containers in parameter list
  plist.AddToList(&tlist);
                  
  // tasks
  MReadTree read("Parameters", inputFile);
  read.DisableAutoScheme();  

  MSrcTranslate srctranslate;
  srctranslate.SetTranslation(xpos,ypos);
  srctranslate.SetRelativeTranslation(kSrcPolicy);

  MSrcRotate    srcrotate;
  srcrotate.SetRAandDEC(fRA,fDEC); 

  MHillasSrcCalc csrc1;   

  MWriteRootFile write(outputFile,"RECREATE");
  write.AddContainer("MHillas"       , "Parameters");
  write.AddContainer("MHillasSrc"    , "Parameters");
  write.AddContainer("MHillasExt"    , "Parameters");
  write.AddContainer("MNewImagePar"  , "Parameters");
  write.AddContainer("MRawEvtHeader" , "Parameters");
  write.AddContainer("MRawRunHeader" , "Parameters");
  write.AddContainer("MConcentration" , "Parameters");
  write.AddContainer("MSrcPosCam"     , "Parameters");

  // include tasks in task list
  tlist.AddToList(&read);
  tlist.AddToList(&srctranslate);
  if(kRotate) 
    tlist.AddToList(&srcrotate);  
  tlist.AddToList(&csrc1);  
  tlist.AddToList(&write);

  // Eventloop
  MEvtLoop evtloop;
  evtloop.SetParList(&plist);
  if (!evtloop.Eventloop(nmaxevents))
    return;  
 
  tlist.PrintStatistics();

}
//-----------------------------------------------------------------------
  
Bool_t readDatacards(TString& filename)
{
  ifstream ifun(filename.Data());
  if(!ifun)
    {
      cout << "File " << filename << " not found" << endl;
      return kFALSE;
    }

  TString word;
  
  while(ifun >> word)
    {
      // skip comments
      if(word[0]=='/' && word[1]=='/')
        {
          while(ifun.get()!='\n'); // skip line
          continue;
        }
      
      // number of events
      if(strcmp(word.Data(),"NEVENTS")==0)
        ifun >> nmaxevents;

      // input file name
      if(strcmp(word.Data(),"INPUTFILES")==0)
        {
          if(inputFile.Length())
            cout << "readDataCards Warning: overriding on-data file name" << endl;
          ifun >> inputFile;
        }

      // output file name
      if(strcmp(word.Data(),"OUTFILE")==0)
        {
          if(outputFile.Length())
            cout << "readDataCards Warning: overriding output file name" << endl;
          ifun >> outputFile;
        }

      // source position
      if(strcmp(word.Data(),"SRCPOS")==0)
        {
          ifun >> xsrcpos;
          ifun >> ysrcpos;        
        }

      // source celestial coordinates 
      if(strcmp(word.Data(),"SRCCOORDS")==0)
        {
          ifun >> fRA;
          ifun >> fDEC;   
        }

      // source movement policy flag
      if(strcmp(word.Data(),"SRCABS")==0)
        ifun >> kSrcPolicy;       

      // rotation flag
      if(strcmp(word.Data(),"ROTFLAG")==0)
        ifun >> kRotate;
    }
  
  // check compulsory values
  if(!inputFile.Length())
    {
      cout << "No on-data file name specified" << endl;
      return kFALSE;
    }

  if(!outputFile.Length())
    {
      cout << "No output file name specified" << endl;
      return kFALSE;
    }

  if(xsrcpos==0 && ysrcpos==0)
    {
      cout << "Source position is center of the camera (as in input file)" << endl;
      return kFALSE;
    }

  // Dump read values
  cout << "************************************************" << endl;
  cout << "* Datacards read from file " << filename << endl;
  cout << "************************************************" << endl;
  cout << "Maximum number of input events: " << nmaxevents << endl;
  cout << "Input file name(s): " << inputFile << endl;
  cout << "Output file name: " << outputFile << endl;
  cout << "Source position (degrees) X=" << xsrcpos << ", Y="<<ysrcpos;
  if(kSrcPolicy)
    cout << " (RELATIVE TO INITIAL SOURCE POSITION)" << endl;
  else
    cout << " (ABSOLUTE POSITION IN THE CAMERA)" << endl;
  cout << "De-rotation flag " << kRotate << endl;
  cout << "Source celestial coordiantes (rad): RA = " << fRA << ", DEC = " << fDEC << endl;
  cout << "***********" << endl << endl;

  return kTRUE;
}