//____________________________________________________________________
// 
// BrChkvPedCalModule
//
// This module is devoted to the Cherenkov detector pedestal calibration 
// It initializes the calibration parameter (BrChkvCalibration)
// Fill histograms during the event loops
// fit them with a gaussian during the finish method and put the
// result in the calibration parameter object
//
// depending on the user's desire, this module can save the stuff into
// an ascii file.
// if you want to commit from this ascii file, it's also possible 
// (cf ChkvCommitCal.C)

//____________________________________________________________________
//
// $Id: BrChkvPedCalModule.cxx,v 1.5 2001/10/23 20:51:23 ouerdane Exp $
// $Author: ouerdane $
// $Date: 2001/10/23 20:51:23 $
// $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov>
//

#ifndef WIN32
#include <iostream>
#include <iomanip>
#include <fstream>
#else
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#endif

#ifndef ROOT_TF1
#include "TF1.h"
#endif
#ifndef ROOT_TString
#include "TString.h"
#endif
#ifndef ROOT_TDirectory
#include "TDirectory.h"
#endif
#ifndef BRAT_BrChkvPedCalModule
#include "BrChkvPedCalModule.h"
#endif
#ifndef BRAT_BrDataTable
#include "BrDataTable.h"
#endif
#ifndef BRAT_BrC1Dig
#include "BrC1Dig.h"
#endif
#ifndef BRAT_BrRichDig
#include "BrRichDig.h"
#endif
#ifndef BRAT_BrParameterDbManager
#include "BrParameterDbManager.h"
#endif
#ifndef BRAT_BrCalibrationManager
#include "BrCalibrationManager.h"
#endif
#ifndef BRAT_BrRunInfoManager
#include "BrRunInfoManager.h"
#endif
#ifndef BRAT_BrRunInfo
#include "BrRunInfo.h"
#endif

//________________________________________________________________________

ClassImp(BrChkvPedCalModule);

//________________________________________________________________________
 BrChkvPedCalModule::BrChkvPedCalModule()
  : BrChkvCalModule()
{
   //
   // default constructor
   //
  SetState(kSetup);
}

//________________________________________________________________________
 BrChkvPedCalModule::BrChkvPedCalModule(Char_t *Name, Char_t *Title) 
  : BrChkvCalModule(Name, Title)
{
   //
   // constructor for named object
   //
  SetState(kSetup);
  SetWidthLimit(); // default is 10 ADC channels
}

//________________________________________________________________________
 BrChkvPedCalModule::~BrChkvPedCalModule()
{
   //
   // destructor
   //
}

//_______________________________________________________________________
 void BrChkvPedCalModule::Init(){
  //
  // Called once per session
  // Get detectorParameters. Inform the parameterElementManeger about
  // database tables to be used (i.e. filled).

  // Job-level initialisation
  SetState(kInit);

  if(DebugLevel() > 1)
    cout << "Entering BrChkvPedCalModule::Init() for " 
	 << GetName() << endl;


  // fCalibration and fParamsChkv set in base class
  BrChkvCalModule::Init();

  // check calibration element access mode

  Int_t nTubes = fParamsChkv->GetNoTubes();
  BrCalibration::EAccessMode mode = BrCalibration::kRead;
  
  // check if we want to load pedestals from ascii file
  if (fLoadAscii) {
    mode = BrCalibration::kTransitional;
    fCalibration->Use("pedestal", mode, nTubes);
    fCalibration->Use("pedestalWidth", mode, nTubes);
  }
  
  // check if we want to save them (tof ascii file or DB)
  else if (fSaveAscii || fCommitAscii) {
    mode = BrCalibration::kWrite;
    fCalibration->Use("pedestal", mode, nTubes);
    fCalibration->Use("pedestalWidth", mode, nTubes);
  }
}

//________________________________________________________________________
 void BrChkvPedCalModule::DefineHistograms()
{
  // -------------------
  // histograms section
  // -------------------

  if (GetState() != kInit) {
    Stop("DefineHistograms", "Must be called after Init"); 
    return;  
  }

  if (fLoadAscii || fCommitAscii) {
    if (Verbose() > 5)
      Warning("DefineHistograms", "No need to declare histos "
	      "since we only load calibration from ascii file"); 
    return;  
  }
    

  TDirectory* saveDir = gDirectory;

  // create a dir for this detector
  fHistDir = gDirectory->mkdir(Form("%s_Pedestals", GetName()));
  fHistDir->mkdir("Tubes");
  
  fHistDir->cd("Tubes");
  // create histos and store them into list
  
  fAdc = new TH1F* [fParamsChkv->GetNoTubes()];
  for(Int_t t = 1; t <= fParamsChkv->GetNoTubes(); t++) {
    fAdc[t-1] = new TH1F(Form("ped%02d", t),
			 Form("%s: Pedestal tube %d", GetName(), t),
			 500, 0, 500);
  }
  
  fHistDir->cd();
  
  // create summaries:
  fSum = new TH1F("peds", "Pedestal summary", fParamsChkv->GetNoTubes(), 
		  0.5, fParamsChkv->GetNoTubes()+.5);
  fSum->SetMarkerStyle(25);
  fSum->SetMarkerColor(2);
  fSum->SetMarkerSize(0.9);

  
  fSumW = new TH1F("pedw", "Pedestal Width summary", 
		   fParamsChkv->GetNoTubes(), 
		   0.5, fParamsChkv->GetNoTubes()+.5);
  fSumW->SetMarkerStyle(25);
  fSumW->SetMarkerColor(2);
  fSumW->SetMarkerSize(0.9);

  // comes back to previous dir
  gDirectory = saveDir;
}

//________________________________________________________________________
 void BrChkvPedCalModule::Begin()
{

  // Run-level initialisation
  SetState(kBegin);

  if (fLoadAscii) {
    ReadAscii();
    return;
  }
  
  // check if histos are booked
  if (!HistBooked())
    if (!fCommitAscii) {
      Abort("Begin", "must book histograms for calibration to work"); 
      return;
    }

  for (Int_t t = 1; t <= fParamsChkv->GetNoTubes(); t++) {
    fCalibration->SetPedestal(t, BrChkvCalibration::kCalException);
    fCalibration->SetPedestalWidth(t, BrChkvCalibration::kCalException);
  }  
}


//________________________________________________________________________
 void BrChkvPedCalModule::Event(BrEventNode* inNode, BrEventNode* outNode)
{
  //
  // Fill the histograms that hold the Pedestal spectra 
  // for syncronization events
  // 
  // Required input e.g. : BrDataTable "DigC1"
  //
  // 
  // Per event method
  
  SetState(kEvent);

  if(DebugLevel() > 1)
    cout << "Entering Event() in BrChkvPedCalModule for " 
	 << GetName() << endl;

  if (fLoadAscii || fCommitAscii)
    return;

  // no base class for BrC1Dig and BrRichDig...
  // have to split code (will fixe that soon)

  if (strcmp(GetName(), "C1") == 0) {
    BrC1Dig* c1 = (BrC1Dig*)inNode->GetObject(Form("Dig%s",GetName()));
    if (!c1) {
      Warning("Event", "No C1 data for this event");
      return;
    }
    
    for(Int_t tube = 1; tube <= fParamsChkv->GetNoTubes(); tube++) 
      fAdc[tube-1]->Fill(c1->GetAdc(tube)); // DO: have to fixe this
    // numbering in BrRawDataInput (like for RICH)
  }
  
  else if (strcmp(GetName(), "RICH") == 0) {
    BrRichDig* rich = (BrRichDig*)inNode->GetObject(Form("Dig%s",GetName()));
    if (!rich) {
      Warning("Event", "No RICH data for this event");
      return;
    }
    
    for(Int_t tube = 1; tube <= fParamsChkv->GetNoTubes(); tube++) 
      fAdc[tube-1]->Fill(rich->GetAdc(tube));
  }
  
  
}

//________________________________________________________________________
 void BrChkvPedCalModule::Finish()
{
  //
  // called once per job 
  //
  
  // Job-level finalisation
  SetState(kFinish);
  
  // if load ascii mode
  if (fLoadAscii) 
    return;

  // if commit mode
  if (fCommitAscii) {
    ReadAscii();
    return;
  }
  
  cout << " Start fitting procedure for " << GetName()
       << " with simple gaussian..."
       << endl;

  if (Verbose())
    cout << " Tube | FitMean | FitSigma | Mean   |  RMS   |"
	 << " DMean  | DSigma " << endl;

  TF1* fit = new TF1("fit", "gaus", 0., 1000.);

  for(Int_t i = 1; i <= fParamsChkv->GetNoTubes() ; i++){
    if (fAdc[i-1]->GetEntries() == 0)
      continue;

    // fit top tube pedestals with gaussian
    
    Float_t mean     = Float_t(fAdc[i-1]->GetMean());
    Float_t pedWidth = Float_t(fAdc[i-1]->GetRMS());
    
    fit->SetParameters(fAdc[i-1]->GetMaximum(), mean, pedWidth);
    fAdc[i-1]->Fit("fit", "Q0", "", 
		   mean - 2*pedWidth, mean + 2*pedWidth);
    
    TF1* f = (TF1*)(fAdc[i-1]->GetListOfFunctions()->Last());
    
    Double_t meanFit  = f->GetParameter(1);    
    Double_t sigmaFit = f->GetParameter(2);
    if (Verbose())
      cout << setw(3)  << i
	   << setw(10) << setprecision(4) << meanFit 
	   << setw(10) << setprecision(4) << sigmaFit
	   << setw(10) << setprecision(4) << mean 
	   << setw(10) << setprecision(4) << pedWidth
	   << setw(10) << setprecision(4) << meanFit - mean
	   << setw(10) << setprecision(4) << sigmaFit - pedWidth 
	   << endl;

    if (TMath::Abs(meanFit - mean) > 20) {
      meanFit  = mean;
      sigmaFit = pedWidth;
      if (Verbose())
	cout << " *** Deviation too large, probably bad fit"
	     << " Take histo MEAN and RMS instead" << endl;
    }

    // set parameters for eventual commit into DB
    
    if (sigmaFit < fWidthLimit) {
      fCalibration->SetPedestal(i, meanFit);
      fCalibration->SetPedestalWidth(i, sigmaFit);
    }

    else 
      if (Verbose())
	cout << " >>>>>> BAD CAL: Width exceeding user's limit <<<<< n";
    
    fSum->Fill(i, meanFit);
    fSumW->Fill(i, sigmaFit);
    
  }
  
  fCalibration->SetComment("pedestal","Generated by BrChkvPedCalModule: "
			   "it with a gaussian function");
  fCalibration->SetComment("pedestalWidth","Generated by BrChkvPedCalModule: "
			   "fit with a gaussian function");
  
  if (fSaveAscii)
    SaveAscii();

}

//____________________________________________________________________
 void BrChkvPedCalModule::SaveAscii() 
{

  // save pedestal to ascii file
    
  BrRunInfoManager* runMan = BrRunInfoManager::Instance();
  Int_t* runs = runMan->GetRunNumbers();
  Int_t  nrun = runMan->GetNumberOfRuns();
  
  BrChkvCalModule::SaveAscii();

  ofstream file(fCalibFile.Data(), ios::out);
  
  file << "****************************************** " << endl;
  file << "*  Calibration for Chkv detector " << GetName() << endl;
  file << "*     Pedestal Calibration " << endl;
  file << "*     Used events from run(s) ";
  for (Int_t r = 0; r < nrun; r++) file << runs[r] << " ";
  file << endl;
  file << "****************************************** " <<endl;
  file << "*" << endl;    
  file << "* tube |  pedestal  | width  " << endl; 
  file << "* ---------------------------" << endl << endl;
  
  for (Int_t tube = 1; tube <= fParamsChkv->GetNoTubes(); tube++)
    file << setw(5) << tube << setw(13) 
	 << fCalibration->GetPedestal(tube) << setw(13)
	 << fCalibration->GetPedestalWidth(tube) << endl;
  
  
  file << "* ------------------------------------" << endl << endl;

  cout << GetName() << " pedestals saved to file " 
       << fCalibFile.Data() << endl << endl;
}


//____________________________________________________________________
 void BrChkvPedCalModule::ReadAscii() 
{

  // save pedestal to ascii file
  BrChkvCalModule::ReadAscii();

  ifstream file(fCalibFile.Data(), ios::in);

  if (!file) {
    Stop("ReadAscii", "File %s was not found", fCalibFile.Data());
    return;
  }

  Float_t ped, wid;
  Int_t tube;
  Char_t comment[256];
  
  file.getline(comment, 256);
  while(comment[0] == '*')
    file.getline(comment, 256);

  if (DebugLevel() > 1) 
    cout << " ---- Pedestals and widths ---- " << endl;
   
  for (Int_t i = 1; i <= fParamsChkv->GetNoTubes(); i++) {
    file >> tube >> ped >> wid;
    if (DebugLevel() > 1) 
      cout << " tube" << setw(3) << tube 
	   << setw(12) << ped << setw(12) << wid << endl;
    
    fCalibration->SetPedestal(tube, ped);
    fCalibration->SetPedestalWidth(tube, wid);
  }
  
  fCalibration->SetComment("pedestal","Generated by BrChkvPedCalModule: "
			   "fit with a gaussian function");
  fCalibration->SetComment("pedestalWidth","Generated by BrChkvPedCalModule: "
			   "fit with a gaussian function");

}

//____________________________________________________________________
 void BrChkvPedCalModule::Print(Option_t* option) const
{
  // Print module information
  // See BrModule::Print for options.
  // In addition this module defines the Option:
  // <fill in here>

  TString opt(option);
  opt.ToLower(); 
  
  BrModule::Print(option); 
  if (opt.Contains("d")) 
   cout << endl 
         << "  Original author: Djamel Ouerdane" << endl
         << "  Last Modifications: " << endl 
         << "    $Author: ouerdane $" << endl  
         << "    $Date: 2001/10/23 20:51:23 $"   << endl 
         << "    $Revision: 1.5 $ " << endl  
         << endl 
         << "-------------------------------------------------" << endl;
}


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

// $Log: BrChkvPedCalModule.cxx,v $
// Revision 1.5  2001/10/23 20:51:23  ouerdane
// Updated modules ala tof or bb, with Set[Save,Commit,Load]Ascii, etc
//
// Revision 1.4  2001/10/08 11:27:33  cholm
// Changed to use new DB access classes
//
// Revision 1.3  2001/08/09 08:58:45  ekman
// Changed tube numbering of C1 - start at 1.
//
// Revision 1.2  2001/08/03 10:11:22  ouerdane
// Started C1 tube numbering from 0
// (note: I still think it should be changed in BrRawDataInput)
// Removed TList member and put histogram members instead
//
// Revision 1.1  2001/07/20 16:05:37  ouerdane
// Added new directory modules/calib/chkv
// Contains:
//   BrChkvCalModule : base class for cherenkov detector calibration modules
//   BrChkvPedCalModule : pedestal calibration module
// stuff for compilation
//   Makefile.am
//   Include.h
//   LinkDef.h
//

Copyright ; 2002 BRAHMS Collaboration <brahmlib@rcf.rhic.bnl.gov>
Last Update on by

Validate HTML
Validate CSS