// $Id: BrZdcVertexCalModule.cxx
// $Author:
// $Date:
// $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov>
//

#if !defined BRAT_BrZdcVertexCalModule
#include "BrZdcVertexCalModule.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 ROOT_TDirectory
#include "TDirectory.h"
#endif

#if !defined ROOT_TF1
#include "TF1.h"
#endif

#ifndef ROOT_TH1
#include "TH1.h"
#endif

#ifndef ROOT_TH2
#include "TH2.h"
#endif

#include <fstream.h>

#if !defined BRAT_BrZdcRdo
#include "BrZdcRdo.h"
#endif

#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 (BrZdcVertexCalModule);

 BrZdcVertexCalModule::BrZdcVertexCalModule ()
{
    SetState (kSetup);

    fVtxOffset= 0;
    fVtxSigma= 0;

    fZdcVsTpc= 0;
    fZdcVsTpc2d= 0;
}



 BrZdcVertexCalModule::BrZdcVertexCalModule (const Char_t *name, const Char_t *title):BrVertexCalModule (name, title)
{
    SetState (kSetup);

    fVtxOffset= 0;
    fVtxSigma= 0;

    fZdcVsTpc= 0;
    fZdcVsTpc2d= 0;

    SetDefaultParameters ();
}



 void BrZdcVertexCalModule::SetDefaultParameters ()
{
    SetFitWindow ();
}



 BrZdcVertexCalModule::~BrZdcVertexCalModule ()
{

}



 void BrZdcVertexCalModule::DefineHistograms ()
{
    if (fCommitAscii || fLoadAscii) return;
    if (GetState ()!= kInit) return;

    TDirectory *saveDir= gDirectory;

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

    fZdcVsTpc= new TH1F ("Zdc", "ZDC vertex- TPM1 track vertices", 600, -300.0, 300.0);
    fZdcVsTpc2d= new TH2F ("Zdc2d", "Zdc vertex- TPM1 track vertices", 200, -100.0, 100.0, 200, -100.0, 100.0);

    gDirectory= saveDir;
}



 void BrZdcVertexCalModule::Init ()
{
    SetState (kInit);

    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 ("zdcVertexOffset", mode, 1);

    if (!fLoadAscii && !fCommitAscii)
    {
        //----- 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 BrZdcVertexCalModule::Begin ()
{
    SetState (kBegin);

    if (fLoadAscii) { ReadAscii (); return; }
    if (fCommitAscii) return;
    if (!HistBooked ()) { Abort ("Begin", "MUST book histograms for calibration to work"); return; }

    BrVertexCalibration::EBrCalException bad= BrVertexCalibration::kCalException;
}



 void BrZdcVertexCalModule::Event (BrEventNode *inNode, BrEventNode *outNode)
{
    SetState(kEvent);

    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;
    }

    //----- get ZDC vertex

    BrZdcRdo *rdoZDC= (BrZdcRdo*) inNode->GetObject ("RdoZDC");
    if (!rdoZDC) return;

    Float_t z= rdoZDC->GetZ ();

    if (TMath::Abs (z) > 200) return;

    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);

        fZdcVsTpc->Fill (z- proj (2));
        fZdcVsTpc2d->Fill (proj (2), z);
    }
}



 void BrZdcVertexCalModule::Finish ()
{
    SetState (kFinish);

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

    TF1 *fit= new TF1 ("fit", "gaus", 0.0, 1000.0);

    //----- fit distribution with gaussian

    if (fZdcVsTpc->GetEntries () > 100)
    {
        Int_t binAtMax= fZdcVsTpc->GetMaximumBin (); // bin number at max
        Axis_t maxZ= fZdcVsTpc->GetBinCenter (binAtMax);

        fit->SetParameters (fZdcVsTpc->GetMaximum (), maxZ, fZdcVsTpc->GetRMS ());

        fZdcVsTpc-> Fit ("fit", "Q0", "", maxZ- fFitWindow, maxZ+ fFitWindow);

        Float_t mean= fit->GetParameter (1);
        Float_t sigma = fit->GetParameter (2);

        fVtxOffset= mean;
        fVtxSigma= sigma;

        BrVertexCalibration::EBrSubDetector sub= BrVertexCalibration::kSubZdc;
        fCalibration->SetVertexOffset (BrVertexCalibration::kZdc, sub, mean, sigma);

        fCalibration-> SetComment ("zdcVertexOffset", "Generated by  BrZdcVertexCalModule: fit with a gaussian function on Z0- TPM1 track Z");

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


    if (Verbose ())
    {
        cout << "ZDC offset:	" << fCalibration->GetZdcOffset () << endl;
        cout << "Sigma:	" << fCalibration->GetZdcSigma () << endl;
    }

    if (fSaveAscii) SaveAscii ();
}



 void BrZdcVertexCalModule::SaveAscii ()
{
    //----- save z-vertex offset and it's sigma to an ASCII file

    Int_t i;
    Float_t par1, par2;

    BrRunInfoManager *runMan= BrRunInfoManager::Instance ();
    const BrRunInfo *run= runMan->GetCurrentRun ();

    if (run->GetRunNo ()== -1) { Abort ("SaveAscii", "RunInfo has run number = -1"); return; }

    BrVertexCalModule::SaveAscii ();
    ofstream file (fCalibFile.Data (), ios::out);

    file.setf (ios::left, ios::adjustfield);
    file.setf (ios::showpoint);
    file.precision (5);
    file.setf (ios::fixed, ios::floatfield);

    file << "***************************************************" << endl;
    file << "*" << endl;
    file << "*" << "	" << GetName () << " z-vertex offset: events from run #" << run->GetRunNo () << endl;
    file << "*" << endl;
    file << "***************************************************" << endl;
    file << "*" << endl;
    file << "*" << "	";

    file.width (8);
    file << "Mod.#";

    file.width (15);
    file << "ZdcOffset";

    file.width (15);
    file << "OffsetSigma" << endl;

    file << "*" << endl;

    par1= fCalibration->GetZdcOffset ();
    par2= fCalibration->GetZdcSigma ();

    file << "	";
    file.width (8);
    file << "1";

    file.width (15);
    file << par1;

    file.width (15);
    file << par2 << endl;
}



 void BrZdcVertexCalModule::ReadAscii ()
{
    BrVertexCalModule::SaveAscii ();
    ifstream file (fCalibFile.Data (), ios::in);

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

    Int_t module;
    Float_t par1, par2;
    Char_t c;

    do
    {
        file.get (c);
        if (c== '*') file.ignore (256, 'n');
    }
    while (c== '*');

    file >> module >> par1 >> par2;

    fCalibration->SetZdcOffset (par1, par2);
    fCalibration->SetComment ("zdcVertexOffset", "ZDC vertex offset relative to TPM1 vertex");
}



 void BrZdcVertexCalModule::Print (Option_t *option)  const
{

}

//----------------------CVS log info
//
// $log$

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

Validate HTML
Validate CSS