source: trunk/MagicSoft/Mars/manalysis/MApplyPadding.cc@ 1683

/* ======================================================================== *\
!
! *
! * 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 <magicsoft@rwagner.de> 10/2002
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MApplyPadding                                                          //
//                                                                         //
//  This task applies padding to a given Sigmabar target value.            //
//  The task checks whether the data stream it is applied to has to be     //
//  padded or not and if so, Gaussian noise with the difference in Sigma   //
//  is produced and added to the particular event. Number of photons and   //
//  error on this number is altered.                                       //
//                                                                         //
//  There are three ways to define the sigmabar value to which the events  //
//  are padded:                                                            //
//                                                                         //
//  1) Set a fixed level with SetTargetLevel                               //
//                                                                         //
//  2) Give a TH1D* which defines the Sigmabar as function of Theta        //
//     with SetDefiningHistogram                                           //
//     The given histogram should have the same binning in Theta as        //
//     is used in the analysis                                             //
//                                                                         //
//  3) Do nothing, then PreProcess will try to read in a (workaround!)     //
//     propriety format ASCII database for the CT1 test.                   //
//     the name of this file is set by SetDatabaseFile                     //
//     Better know what you are doing or use methods 1 or 2.               //
//                                                                         //
//  This implementation is still PRELIMINARY and requires some workarounds //
//  put in SPECIFICALLY FOR THE CT1 TESTS, since a database to access is   //
//  missing. It is not the FINAL MAGIC VERSION.                            //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MApplyPadding.h"

#include <stdio.h>

#include "TH1.h"
#include "TH2.h"
#include "TRandom.h"

#include "MBinning.h"
#include "MSigmabar.h"
#include "MMcEvt.hxx"
#include "MLog.h"
#include "MLogManip.h"
#include "MParList.h"
#include "MGeomCam.h"
#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"
#include "MPedestalPix.h"

ClassImp(MApplyPadding);

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MApplyPadding::MApplyPadding(const char *name, const char *title) : fRunType(0), fGroup(0), fUseHistogram(kTRUE), fFixedSigmabar(0.0)
{
  fName  = name  ? name  : "MApplyPadding";
  fTitle = title ? title : "Task to apply padding";
  Print();
}

// --------------------------------------------------------------------------
//
// Destructor. 
//
MApplyPadding::~MApplyPadding()
{
  //nothing yet
}

// --------------------------------------------------------------------------
//
// You can provide a TH1D* histogram containing the target Sigmabar in
// bins of theta. Be sure to use the same binning as for the analysis
//
Bool_t MApplyPadding::SetDefiningHistogram(TH1D *histo)
{
  fHSigmabarMax = histo;
  return kTRUE;
}


// --------------------------------------------------------------------------
//
//  check if MEvtHeader exists in the Parameter list already.
//  if not create one and add them to the list
//
Bool_t MApplyPadding::PreProcess(MParList *pList)
{
  fRnd = new TRandom3(0);

  fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
  if (!fMcEvt)
    {
       *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
       return kFALSE;
     }
  
   fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
   if (!fPed)
     {
       *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
       return kFALSE;
     }
  
   fCam = (MGeomCam*)pList->FindObject("MGeomCam");
   if (!fCam)
     {
       *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
       return kFALSE;
     }
  
   fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
   if (!fEvt)
     {
       *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
       return kFALSE;
     }

   fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
   if (!fSigmabar)
     {
       *fLog << dbginf << "MSigmabar not found... aborting." << endl;
       return kFALSE;
     }
   
   // Get Theta Binning  
   MBinning* binstheta  = (MBinning*)pList->FindObject("BinningTheta");
   if (!binstheta)
     {
       *fLog << err << dbginf << "BinningTheta not found... aborting." << endl;
       return kFALSE;      
     }

   // Create fSigmabarMax histogram
   // (only if no fixed Sigmabar target value or a histogram have already been 
   // provided)
   if ((!fUseHistogram) && (fHSigmabarMax==NULL)) {
     
     fHSigmabarMax = new TH1D();
     fHSigmabarMax->SetNameTitle("fHSigmabarMax","Sigmabarmax for this analysis");
     TAxis &x = *fHSigmabarMax->GetXaxis();
#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
     TString xtitle = x.GetTitle();
#endif
     fHSigmabarMax->SetBins(binstheta->GetNumBins(), 0, 1);
     // Set the binning of the current histogram to the binning
     // in one of the two given histograms
     x.Set(binstheta->GetNumBins(), binstheta->GetEdges());
#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
     x.SetTitle(xtitle);
#endif
     
     // -------------------------------------------------
     // read in SigmabarParams
     // workaround--proprietary file format--CT1test only BEGIN
     // -------------------------------------------------
     
     FILE *f;
     if( !(f =fopen(fDatabaseFilename, "r")) ) {
       *fLog << err << dbginf << "Database file " << fDatabaseFilename << "was not found... (specify with MApplyPadding::SetDatabaseFile) aborting." << endl;
       return kFALSE;  
     }
     char line[80];    
     Float_t sigmabarMin, sigmabarMax, thetaMin, thetaMax, ra, dec2;
     Int_t type, group, mjd, nr;
     while ( fgets(line, sizeof(line), f) != NULL) {
       if ((line[0]!='#')) {    
         sscanf(line,"%d %d %f %f %d %d %f %f %f %f",&type, &group, &ra, &dec2, &mjd, &nr, &sigmabarMin,&sigmabarMax,&thetaMin,&thetaMax);
         if ((group==fGroup)||(type==1)) //selected ON group or OFF         
           {    
             // find out which bin(s) we have to look at
             for (Int_t i=fHSigmabarMax->GetXaxis()->FindBin(thetaMin); 
                  i<fHSigmabarMax->GetXaxis()->FindBin(thetaMax)+1; i++) 
               if (sigmabarMax > fHSigmabarMax->GetBinContent(i)) 
                 fHSigmabarMax->SetBinContent(i, sigmabarMax);    
           }
       }
     }//while
  
     // workaround--proprietary file format--CT1test only END

     //   fHSigmabarMax->DrawClone();     
     //   fTest = new TH2D("fTest", "Test if padding works", 201, -0.05, 2.05, 201, -0.05, 2.05);

   } //!fUseHistogram
   return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculate Sigmabar for current event
// Then apply padding
// 
// 1) have current event
// 2) get sigmabar(theta)
// 3) pad event
//
Bool_t MApplyPadding::Process()
{
  // Calculate sigmabar of event
  fSigmabar->Calc(*fCam, *fPed);
  Double_t mySig = fSigmabar->GetSigmabar();

  // Get sigmabar which we have to pad to
  Double_t otherSig;
  if (fUseHistogram) {
    Int_t binNumber = fHSigmabarMax->GetXaxis()->FindBin(fMcEvt->GetTheta()*kRad2Deg);
    otherSig = fHSigmabarMax->GetBinContent(binNumber);
  } else {
    otherSig = fFixedSigmabar;
  }

  // Determine quadratic difference other-mine
  Double_t quadraticDiff = otherSig*otherSig - mySig*mySig;

  if (quadraticDiff < 0) {
    *fLog << err << dbginf << "Event has higher Sigmabar="<<mySig<<" than Sigmabarmax="<<otherSig << " ...Skipping this event" <<endl;
    return kCONTINUE; //skip
  }

  if (quadraticDiff == 0) return kTRUE; //no padding necessary.
  
  // Pad if quadratic difference > 0 
  if (quadraticDiff > 0) {    

   MPedestalCam newPed;
   newPed.InitSize(fPed->GetSize());

   const UInt_t npix = fPed->GetSize(); // Total number of pixels
   for (UInt_t i=0; i<npix; i++) {
     MCerPhotPix pix = fEvt->operator[](i);
     if (!pix.IsPixelUsed())
       continue;
     pix.SetNumPhotons(pix.GetNumPhotons() +
                       sqrt(quadraticDiff)*
                       fRnd->Gaus(0.0, 1.0)/
                       fCam->GetPixRatio(pix.GetPixId())
                       );
     // error: add sigma of padded noise quadratically
     Double_t error = pix.GetErrorPhot();
     pix.SetErrorPhot(sqrt(error*error + quadraticDiff));
     
     MPedestalPix ppix = fPed->operator[](i);
     MPedestalPix npix;
     npix.SetSigma(sqrt(ppix.GetSigma()*ppix.GetSigma() + quadraticDiff));
     newPed[i]=npix;
    } //for
   // Calculate Sigmabar again and crosscheck
   fSigmabar->Calc(*fCam, newPed);
   //mySig = fSigmabar->GetSigmabar();
   // fTest->Fill(otherSig,mySig);
   return kTRUE;
  } //if 
  return kFALSE;
}

Bool_t MApplyPadding::PostProcess()
{
  //  fTest->DrawClone();
  return kTRUE;
}