//____________________________________________________________________
// 
// Module for pedestal calibration of beam-beam counters 
// 

//____________________________________________________________________
//
// $Id: BrInelVertexCalModule.cxx,v 1.3 2002/08/30 16:19:54 hagel Exp $
// $Author: hagel $
// $Date: 2002/08/30 16:19:54 $
// $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov>
//

#if !defined BRAT_BrInelVertexCalModule
#include "BrInelVertexCalModule.h"
#endif
#if !defined BRAT_BrEventNode
#include "BrEventNode.h"
#endif
#if !defined BRAT_BrVertexCalibration
#include "BrVertexCalibration.h"
#endif
#if !defined ROOT_TString
#include "TString.h"
#endif
#if !defined BRAT_BrDetectorVolume
#include "BrDetectorVolume.h"
#endif
#if !defined BRAT_BrDataTable
#include "BrDataTable.h"
#endif
#if !defined BRAT_BrInelVertex
#include "BrInelVertex.h"
#endif

#if !defined ROOT_TDirectory
#include "TDirectory.h"
#endif
#if !defined ROOT_TF1
#include "TF1.h"
#endif
#if !defined ROOT_TString
#include "TString.h"
#endif
#ifndef ROOT_TH1
#include "TH1.h"
#endif
#ifndef ROOT_TH2
#include "TH2.h"
#endif

#include <fstream.h>
#if !defined BRAT_BrIostream
#include "BrIostream.h"
#endif

#if !defined BRAT_BrRunInfoManager
#include "BrRunInfoManager.h"
#endif
#if !defined BRAT_BrDetectorTrack
#include "BrDetectorTrack.h"
#endif
#if !defined BRAT_BrPlane3D
#include "BrPlane3D.h"
#endif
#if !defined BRAT_BrVector3D
#include "BrVector3D.h"
#endif
#if !defined BRAT_BrLine3D
#include "BrLine3D.h"
#endif
#if !defined BRAT_BrGeometryDbManager
#include "BrGeometryDbManager.h"
#endif

//____________________________________________________________________
ClassImp(BrInelVertexCalModule);

//____________________________________________________________________
 BrInelVertexCalModule::BrInelVertexCalModule()
{
  // Default constructor. DO NOT USE
  SetState(kSetup);

  fVtxOffset = 0;
  fInelVsTpc   = 0;
  fInelVsTpc2d = 0;

  fVertexModule = 0;
  fVtxOffset = 0;
  fVtxSigma  = 0;
}

//____________________________________________________________________
 BrInelVertexCalModule::BrInelVertexCalModule(const Char_t* name, const Char_t* title)
  : BrVertexCalModule(name, title)
{
  // Named Constructor
  SetState(kSetup);

  fVtxOffset = 0;
  fInelVsTpc   = 0;
  fInelVsTpc2d = 0;

  fVertexModule = new BrInelVertexModule( "Inel", "Inel Vertex Module");

  SetDefaultParameters();

}

//____________________________________________________________________
 void BrInelVertexCalModule::SetDefaultParameters()
{
  // set default parameters

  fVertexModule->SetCorrectTdc(kTRUE);
  fVertexModule->SetOffsetCalib(kTRUE);

  SetFitWindow();       // default is 16 cm (2 sigmas) 
}

//____________________________________________________________________
 BrInelVertexCalModule::~BrInelVertexCalModule()
{
  // destructor
  if (fVertexModule)
    delete fVertexModule;
}

//____________________________________________________________________
 void BrInelVertexCalModule::DefineHistograms()
{
  // Define histograms. They are:
  // <fill in here>

  if (fCommitAscii || fLoadAscii) {
    if (Verbose() > 10)
      Warning("DefineHistograms", "No need to be here if you "
	      "load or commit an ascii calibration");
    return;
  }

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

  TDirectory* saveDir = gDirectory; 

  // book sub module histos

  fVertexModule->Book();

  fHistDir = gDirectory->mkdir("Inel_VertexCal"); 
  fHistDir->cd();

  // Make histograms here 
  fInelVsTpc = 
    new TH1F("Inel", "inel vertex - TPM1 track vertices", 
	     200, -100, 100);

  fInelVsTpc2d = 
    new TH2F("bInel2d", "Inel vertex - TPM1 track vertices", 
	     50, -100, 100, 100, -100, 100);
  
  gDirectory = saveDir;
}

//____________________________________________________________________
 void BrInelVertexCalModule::Init()
{
  // Job-level initialisation
  SetState(kInit);

  //
  // Called once per session
  // Get detectorParameters. Inform the parameterElementManeger about
  // database tables to be used (i.e. filled).

  if(DebugLevel() > 2)
    cout << "Entering BrInelVertexCalModule::Init() for " 
	 << GetName() << endl;

  // initializing base class

  BrVertexCalModule::Init();
  BrCalibration::EAccessMode mode = BrCalibration::kRead;
  
  if(fLoadAscii)
    mode = BrCalibration::kTransitional;
  else if (fSaveAscii || fCommitAscii)
    mode = BrCalibration::kWrite;
  else {
    Abort("Init", "You MUST either save to ascii, load from ascii "
	  "or commit from ascii. Use the proper setter! " );
    return;
  }

  fCalibration->Use("inelVertexOffset", mode, 2);

  if (!fLoadAscii && ! fCommitAscii) {
    // initialize hit and vertex modules
    fVertexModule->SetOffsetCalib(kTRUE);
    fVertexModule->Init();
    fVertexModule->SetVerbose(Verbose());
    fVertexModule->SetDebugLevel(DebugLevel());
    
    // get TPM1 volume
    BrGeometryDbManager* geo = BrGeometryDbManager::Instance();
    fTpm1Volume = (BrDetectorVolume*)geo->
      GetDetectorVolume("BrDetectorVolume", "TPM1");
    if(!fTpm1Volume){
      cout << "No TPM1 volume" << endl;
      Abort("Init","No TPM1 detector");
    }

  }
}

//____________________________________________________________________
 void BrInelVertexCalModule::Begin()
{
  // Run-level initialisation
  SetState(kBegin);
  
  if(fLoadAscii){
    ReadAscii();    
    return;
  }
  
  if (fCommitAscii)
    return;
  
  // check if histos are booked
  if (!HistBooked()) {
    Abort("Begin", "MUST book histograms for calibration to work"); 
    return;
  }

  BrVertexCalibration::EBrCalException bad = 
    BrVertexCalibration::kCalException;
  fCalibration->SetInelOffset(bad, bad);

  fVertexModule->Begin();
  
}

//____________________________________________________________________
 void BrInelVertexCalModule::Event(BrEventNode* inNode, BrEventNode* outNode)
{
  // Per event method
  SetState(kEvent);

  // -------------------------------------------------------------------
  // Get TPM1 tracks and Inel vertex
  // project tracks to longitudinal plane (X = 0 plane)
  // fill histos with Z0 vs track proj Z and difference between those
  // quantities.
  // -------------------------------------------------------------------
  if(DebugLevel() > 1)
    cout << "Entering Event() in BrInelVertexCalModule for " 
	 << GetName() << endl;

  if (fCommitAscii || fLoadAscii)
    return;

  
  // Get TPM1 tracks
  BrDataTable* trks = inNode->GetDataTable("DetectorTrack TPM1");
  if (!trks) {
    if (DebugLevel() > 2) Warning("Event", "No TPM1 tracks for this event");
    return;
  }
  
  // call hit and vertex modules event methods

  fVertexModule->Event(inNode, outNode);
  
  // Get Inel Vertex
  BrInelVertex* bbs = (BrInelVertex*)outNode->GetObject("InelVertex");
  //printf("bbs = %dn",Int_t(bbs));
  if (bbs)
    Process(bbs, trks); 
}

//____________________________________________________________________
 void BrInelVertexCalModule::Process(BrInelVertex* bb, BrDataTable* trks)
{
  // private method
  // loop over tracks and project them to longitudinal plane
  // fill histograms InelZ0 - Z track proj.
  
  // plane at X = 0
  BrPlane3D plane(0,0,0, 0,1,0, 0,0,1);
  for (Int_t t = 0; t < trks->GetEntries(); t++) {
    BrDetectorTrack* trk = (BrDetectorTrack*)trks->At(t);
    BrLine3D trkLine = fTpm1Volume->LocalToGlobal(trk->GetTrackLine());
    BrVector3D proj = plane.GetIntersectionWithLine(trkLine);
    
    fInelVsTpc->Fill(bb->GetZ() - proj(2));
    fInelVsTpc2d->Fill(proj(2), bb->GetZ());
    Info(50,"Process","Inel vertex z = %f, tpm1 vertex z = %fn",bb->GetZ(),proj(2));
  }
}

//____________________________________________________________________
 void BrInelVertexCalModule::Finish()
{
  // Job-level finalisation.
  // The requirements for a succesfull fit at end is
  // a) More than 100 counts in spectrum -- even this is marginal
  // b) If < 2000 counts fit with maximum likelyhood to improve quality
  //    of fit.
  // It is concieveable that 
  //
  SetState(kFinish);
  
  if(fLoadAscii)
    return;

  if (fCommitAscii) {
    ReadAscii();
    return;
  }

  fVertexModule->Finish();
  
  TF1* fit = new TF1("fit", "gaus", 0., 1000.);

    // fit distribution with gaussian
  if (fInelVsTpc->GetEntries() > 100){

    
    Int_t binAtMax  = fInelVsTpc->GetMaximumBin(); // bin number at max
    Axis_t maxZ = fInelVsTpc->GetBinCenter(binAtMax);
    
    fit->SetParameters(fInelVsTpc->GetMaximum(), maxZ, fInelVsTpc->GetRMS());
    
    if (fInelVsTpc->GetEntries() > 2000)
      fInelVsTpc->Fit("fit", "Q0", "", maxZ - fFitWindow, maxZ + fFitWindow);
    else
      fInelVsTpc->Fit("fit", "Q0L", "", maxZ - fFitWindow, maxZ + fFitWindow);

    
    Float_t mean  = fit->GetParameter(1);    
    Float_t sigma = fit->GetParameter(2);
    
    fVtxOffset = mean;
    fVtxSigma  = sigma;
    
    BrVertexCalibration::EBrSubDetector sub = BrVertexCalibration::kSubInel;
    fCalibration->SetVertexOffset(BrVertexCalibration::kInel, sub, mean, sigma);
    
    fCalibration->
      SetComment("inelVertexOffset",
		 "Generated by  BrInelVertexCalModule: "
		 "fit with a gaussian function on Z0 - TPM1 track Z");

  }
  else
    Warning("Finish","Fit did not succeed - too few counts");

  if (Verbose())
    cout << " --------- Inel vertices offsets with TPM1 ---------- " << endl
	 << "  Inel offset    : " 
	 << setw(8) << fCalibration->GetInelOffset() 
	 << " +/- " 
	 << setw(7) << fCalibration->GetInelSigma() << " cm " << endl;
  
  if (fSaveAscii)
    SaveAscii();
  
}

//____________________________________________________________________
 void BrInelVertexCalModule::SaveAscii() 
{

  // save pedestal to ascii file

  BrRunInfoManager* runMan = BrRunInfoManager::Instance();
  Int_t* runs = runMan->GetRunNumbers();
  Int_t  nrun = runMan->GetNumberOfRuns(); 
  
  BrVertexCalModule::SaveAscii();

  ofstream file(fCalibFile.Data(), ios::out);
  
  file << "****************************************** " << endl;
  file << "*  Calibration for Inel counter " << GetName() << endl;
  file << "*     Inel vertex offsets along the beam line" << 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 << "*    Inel PP counters " << endl;
  file << "* Offset  -  Sigma  " << endl;
  file << "* --------------------------------------------------------" 
       << endl << endl;
  
  file << setw(10) << setprecision(5) << fCalibration->GetInelOffset()   << "  " 
       << setw(10) << setprecision(5) << fCalibration->GetInelSigma()    << endl;
  file << "Comment:" << endl;
  //  file << fCalibration->GetComment() << endl;
  file << "* --------------------------------------------------------" << endl << endl;
}

//____________________________________________________________________
 void BrInelVertexCalModule::ReadAscii() 
{
  
  // save pedestal to ascii file
  BrVertexCalModule::SaveAscii();

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

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

  Float_t big, small, fast, sbig, ssmall, sfast;
  Char_t comment[256];
  
  file.getline(comment, 256);
  while(comment[0] == '*')
    file.getline(comment, 256);
  
  file >> big >> sbig ;

  if (DebugLevel() > 2) 
    cout << " Values read in the file " << fCalibFile.Data() << endl
	 << "    " << setw(8) << big   << " +/- " << setw(8) << sbig   << endl;

  fCalibration->SetInelOffset(big, sbig);
  fCalibration->SetComment("inelVertexOffset", fComment.Data());
}


//____________________________________________________________________
 void BrInelVertexCalModule::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: Flemming Videbaek" << endl
         << "  Last Modifications: " << endl 
         << "    $Author: hagel $" << endl  
         << "    $Date: 2002/08/30 16:19:54 $"   << endl 
         << "    $Revision: 1.3 $ " << endl  
         << endl 
         << "-------------------------------------------------" << endl;
}

//____________________________________________________________________
//
// $Log: BrInelVertexCalModule.cxx,v $
// Revision 1.3  2002/08/30 16:19:54  hagel
// Changed limits on 2d diagnostic histogram
//
// Revision 1.2  2002/08/15 17:20:32  videbaek
// Increase range in 2D histogram
//
// Revision 1.1  2002/06/07 16:01:53  videbaek
// Include the calibration module for the Inelastic vertex in pp.
//

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

Validate HTML
Validate CSS