source: trunk/MagicSoft/Mars/manalysis/MPadding.cc@ 1885

/* ======================================================================== *\
!
! *
! * 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   <mailto:magicsoft@rwagner.de> 10/2002
!              Wolfgang Wittek <mailto:wittek@mppmu.mpg.de>  01/2003
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MPadding                                                               //
//  (developped from 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. The number of photons,  // 
//  its error and the pedestal sigmas are altered.                         //
//                                                                         //
//  The padding has to be done before the image cleaning because the       //
//  image cleaning depends on the pedestal sigmas.                         //
//                                                                         //
//  There are several ways of defining the sigmabar value to which the     // 
//  events are padded:                                                     //
//                                                                         //
//  1) Set a fixed level (fFixedSigmabar) by calling 'SetTargetLevel'.     //
//                                                                         //
//  2) By calling 'SetDefiningHistogram', give a TH1D histogram            //
//     (fHSigmabarMax) which defines the Sigmabar as a function of Theta.  //
//                                                                         //
//  3) By calling 'SetSigmaThetaHist', give a TH2D histogram               //
//     (fHSigmaTheta) which contains the Sigmabar distribution for the     //
//     different bins in Theta. For a given event, the sigmabar value to   //
//     be used for the padding is thrown from this distribution.           //
//                                                                         //
//  Workaround :                                                           //  
//  If none of these options is specified then PreProcess will try to read // 
//  in a propriety format ASCII database for the CT1 test. The name of     // 
//  this file is set by 'SetDatabaseFile'. From the data in this file a    //
//  TH1D histogram (fHSigmabarMax) is generated.                           //
//                                                                         //
//  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 "MPadding.h"

#include <stdio.h>

#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TProfile.h"
#include "TRandom.h"
#include "TCanvas.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(MPadding);

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MPadding::MPadding(const char *name, const char *title) : fRunType(0), fGroup(0), fFixedSigmabar(0.0)
{
  fName  = name  ? name  : "MPadding";
  fTitle = title ? title : "Task for the padding";

  fFixedSigmabar = 0.0;
  fHSigmabarMax  = NULL;
  fHSigmaTheta   = NULL;

  Print();
}

// --------------------------------------------------------------------------
//
// Destructor. 
//
MPadding::~MPadding()
{
  //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 MPadding::SetDefiningHistogram(TH1D *histo)
{
  fHSigmabarMax  = histo;

  fFixedSigmabar = 0.0;
  fHSigmaTheta   = NULL;

  *fLog << "MPadding::SetDefiningHistogram" << endl;

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// You can provide a TH2D* histogram containing 
//     - the Sigmabar distribution in bins of theta. 
//
Bool_t MPadding::SetSigmaThetaHist(TH2D *histo)
{
  fHSigmaTheta   = histo;

  fFixedSigmabar = 0.0;
  fHSigmabarMax  = NULL;

  *fLog << "MPadding::SetSigmaThetaHist" << endl;

  return kTRUE;
}

// --------------------------------------------------------------------------
//
//
void MPadding::SetTargetLevel(Double_t sigmabar)
{
  fFixedSigmabar = sigmabar;

  fHSigmabarMax  = NULL;
  fHSigmaTheta   = NULL;

  *fLog << "MPadding::SetTargetLevel; use fixed sigmabar : fFixedSigmabar = " 
        << fFixedSigmabar << endl;
}

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

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

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

   // Get Sigma Binning  
   MBinning* binssigma  = (MBinning*)pList->FindObject("BinningSigmabar");
   if (!binssigma)
     {
       *fLog << err << dbginf << "MPadding::PreProcess : BinningSigmabar not found... aborting." << endl;
       return kFALSE;      
     }

   //--------------------------------------------------------------------
   // plot pedestal sigmas for testing purposes
   fHSigmaPedestal = new TH2D("SigPed","padded vs orig. sigma", 
                     100, 0.0, 5.0, 100, 0.0, 5.0);
   fHSigmaPedestal->SetXTitle("orig. Pedestal sigma");
   fHSigmaPedestal->SetYTitle("padded Pedestal sigma");

   // plot no.of photons (before vs. after padding) for testing purposes
   fHPhotons = new TH2D("Photons","Photons after vs.before padding", 
                        100, -10.0, 90.0, 100, -10, 90);
   fHPhotons->SetXTitle("no.of photons before padding");
   fHPhotons->SetYTitle("no.of photons after padding");

   // plot of added NSB
   fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0);
   fHNSB->SetXTitle("no.of added NSB photons");
   fHNSB->SetYTitle("no.of pixels");

   fHSigmaOld = new TH2D();
   fHSigmaOld->SetNameTitle("fHSigmaOld","Sigma before padding");
   MH::SetBinning(fHSigmaOld, binstheta, binssigma);
   fHSigmaOld->SetXTitle("Theta");
   fHSigmaOld->SetYTitle("Sigma");

   fHSigmaNew = new TH2D();
   fHSigmaNew->SetNameTitle("fHSigmaNew","Sigma after padding");
   MH::SetBinning(fHSigmaNew, binstheta, binssigma);
   fHSigmaNew->SetXTitle("Theta");
   fHSigmaNew->SetYTitle("Sigma");


   //************************************************************************
   // Create fSigmabarMax histogram
   // (only if no fixed Sigmabar target value and no histogram have been 
   // provided)
   //
   if ( (fFixedSigmabar==0.0)   && (fHSigmabarMax==NULL) 
        && (fHSigmaTheta==NULL) ) 
   {
        *fLog << "MPadding::PreProcess() : fFixedSigmabar, fHSigmabarMax = " 
          <<  fFixedSigmabar << ",  " << fHSigmabarMax << endl;
     *fLog << "          create fSigmabarMax histogram" << endl;

     fHSigmabarMax = new TH1D();
     fHSigmabarMax->SetNameTitle("fHSigmabarMax","Sigmabarmax for this analysis");
     TAxis &x = *fHSigmabarMax->GetXaxis();
     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());
     x.SetTitle("Theta");
     TAxis &y = *fHSigmabarMax->GetYaxis();     
     y.SetTitle("SigmabarMax");

     // -------------------------------------------------
     // 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 MPadding::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
  
   } //fFixedSigmabar
   //************************************************************************

   return kTRUE;
}

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

  //$$$$$$$$$$$$$$$$$$$$$$$$$$
  mySig = 0.0;
  //$$$$$$$$$$$$$$$$$$$$$$$$$$


  // Get sigmabar which we have to pad to
  //
  Double_t otherSig;
  Double_t Theta  = kRad2Deg*fMcEvt->GetTelescopeTheta();

  // *fLog << "Theta = " << Theta << endl;

  //-------------------------------------------
  // get target sigma for the current Theta from the histogram fHSigmabarMax
  //

  if (fHSigmabarMax != NULL) 
  {
    Int_t binNumber = fHSigmabarMax->GetXaxis()->FindBin(Theta);
    if (binNumber < 1  || binNumber > fHSigmabarMax->GetNbinsX())
    {
      *fLog << err << dbginf 
            << "Theta of current event is beyond the limits, Theta = "  
            << kRad2Deg*fMcEvt->GetTelescopeTheta()
            << "   ...Skipping this event" <<endl;
      return kCONTINUE;
    }
    else
    {
      otherSig = fHSigmabarMax->GetBinContent(binNumber);

      //*fLog << "Theta, binNumber, otherSig = " 
      //      << kRad2Deg*fMcEvt->GetTelescopeTheta() << ",  "
      //      << binNumber << ",  " << otherSig << endl; 
    }
  } 

  //-------------------------------------------
  // for the current Theta, 
  // generate a sigma according to the histogram fHSigmaTheta
  //
  else if (fHSigmaTheta != NULL) 
  {
    Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(Theta);

    if ( binNumber < 1  ||  binNumber > fHSigmaTheta->GetNbinsX() )
    {
      //      *fLog << "MPadding::Process(); binNumber out of range, binNumber = "
      //      << binNumber << ",  Skip event " << endl;
      return kCONTINUE;
    }
    else
    {
      TH1D* fHSigma = 
            fHSigmaTheta->ProjectionY("", binNumber, binNumber, "");
      if ( fHSigma->GetEntries() == 0.0 )
      {
        //*fLog << "MPadding::Process(); projection for Theta bin " << binNumber
        //      << " has no entries,  Skip event" << endl;
        return kCONTINUE;
      }
      else
      {
        otherSig = fHSigma->GetRandom(); 

        //*fLog << "Theta, bin number = " << Theta << ",  " << binNumber 
        //      << ",  otherSig = " << otherSig << endl;
      }
      delete fHSigma;
    }
  } 

  //-------------------------------------------
  // use a fixed target sigma
  //
  else if (fFixedSigmabar != 0.0) 
  {
    otherSig = fFixedSigmabar;
  }

  //-------------------------------------------
  else
  {
    *fLog << "MPadding::Process(); sigmabar for padding not defined" << endl;
    return kFALSE;
  }

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

  //*fLog << "MPadding::Process(); mySig, otherSig = " << mySig << ",  "
  //      << otherSig << endl;

  // Skip event if target sigma is zero
  if (otherSig == 0.0)
  {
    return kCONTINUE;     
  }

  // 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 << "; Theta =" 
          << kRad2Deg*fMcEvt->GetTelescopeTheta()
          << "   ...Skipping this event" <<endl;
    return kCONTINUE; //skip
  }

  if (quadraticDiff == 0) return kTRUE; //no padding necessary.
  

  //-------------------------------------------
  // quadratic difference is > 0, do the padding; 
  //
  Padding(quadraticDiff);

  // Calculate Sigmabar again and crosscheck
  mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);

  //fSigmabar->Print("");

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Do the padding  (mySig ==> otherSig)
// 
Bool_t MPadding::Padding(Double_t quadraticDiff)
{
   const UInt_t npix = fEvt->GetNumPixels();

   // pad only pixels   - which are used (before image cleaning)
   //                   - and for which the no.of photons is != 0.0
   //
   for (UInt_t i=0; i<npix; i++) 
   {   
     MCerPhotPix &pix = fEvt->operator[](i);      
     if ( !pix.IsPixelUsed() )
       continue;

     if ( pix.GetNumPhotons() == 0.0)
     {
       *fLog << "MPadding::Process(); no.of photons is 0 for used pixel" 
             << endl;
       continue;
     }

     Int_t j = pix.GetPixId();
     Double_t Area = fCam->GetPixRatio(j);

     // add additional NSB to the number of photons
     Double_t oldphotons = pix.GetNumPhotons();
     Double_t NSB = sqrt(quadraticDiff*Area) * fRnd->Gaus(0.0, 1.0);
     Double_t newphotons = oldphotons + NSB;
     pix.SetNumPhotons( newphotons );

     fHNSB->Fill( NSB/sqrt(Area) );
     fHPhotons->Fill( newphotons/sqrt(Area), oldphotons/sqrt(Area) );


     // error: add sigma of padded noise quadratically
     Double_t olderror = pix.GetErrorPhot();
     Double_t newerror = sqrt(  olderror*olderror + quadraticDiff*Area );
     pix.SetErrorPhot( newerror );


     MPedestalPix &ppix = fPed->operator[](j);

     //$$$$$$$$$$$$$$$$$$$$$$$$$$
     ppix.SetMeanRms(0.0);
     //$$$$$$$$$$$$$$$$$$$$$$$$$$

     Double_t oldsigma = ppix.GetMeanRms();
     Double_t newsigma = sqrt(  oldsigma*oldsigma + quadraticDiff*Area );
     ppix.SetMeanRms( newsigma );

     fHSigmaPedestal->Fill( oldsigma, newsigma );
     fHSigmaOld->Fill( kRad2Deg*fMcEvt->GetTelescopeTheta(), oldsigma );
     fHSigmaNew->Fill( kRad2Deg*fMcEvt->GetTelescopeTheta(), newsigma );
   } //for

   return kTRUE;
}

// --------------------------------------------------------------------------
//
//
Bool_t MPadding::PostProcess()
{
    TCanvas &c = *(MH::MakeDefCanvas("Padding", "", 600, 900)); 
    c.Divide(2,3);
    gROOT->SetSelectedPad(NULL);


    if (fHSigmabarMax != NULL)
    {
      c.cd(1);
      fHSigmabarMax->DrawClone();     
      fHSigmabarMax->SetBit(kCanDelete);     
    }
    else if (fHSigmaTheta != NULL)
    {
      c.cd(1);
      fHSigmaTheta->DrawClone();     
      fHSigmaTheta->SetBit(kCanDelete);     
    }

    c.cd(3);
    fHSigmaPedestal->DrawClone();
    fHSigmaPedestal->SetBit(kCanDelete);    

    c.cd(5);
    fHPhotons->DrawClone();
    fHPhotons->SetBit(kCanDelete);    

    c.cd(2);
    fHNSB->DrawClone();
    fHNSB->SetBit(kCanDelete);    

    c.cd(4);
    fHSigmaOld->DrawClone();     
    fHSigmaOld->SetBit(kCanDelete);     

    c.cd(6);
    fHSigmaNew->DrawClone();     
    fHSigmaNew->SetBit(kCanDelete);     


  return kTRUE;
}

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