BRAT 2.4.5
Class index
Full class index
brahmlib
BRAHMS
ROOT page 
//____________________________________________________________________
// 
// Class BrTpcTimeCalModule
//
// This module atempts to correct for the row time offsets seen in the
// Tpcs. 
//
// INPUT : The module needs tracks for the TPC, but these _HAVE_ to be
// done using a loose cut on the on the y search width.
// Example :
//    Char_t *tpcNames[4] = {"TPM1", "TPM2", "T1", "T2"};
//
//    for(Int_t i = 1; i < 4; i++) {
//    
//      BrTpcTrackPackage* tpcTrackPackage = 
//        new BrTpcTrackPackage(tpcNames[i], "Hit package");
//
//      mainModule->AddModule(tpcTrackPackage);
//
//      BrTpcTrackFollowModule* trackModule = 
//        (BrTpcTrackFollowModule*)tpcTrackPackage->GetTpcTrackModule();
//      trackModule->SetMaxRowsMissed(2);
//      trackModule->SetSearchWidthX(0.2);
//      trackModule->SetSearchWidthY(10.0);  <-- NB!
//    }
//
// IDEA : The idea is to assume that some rows, in the middle, are
// only slightly affected. We use those to refit the track and
// calculate the deviation in the rows close to each end as well as in
// the middle.
// 
// HOW TO USE MODULE : Example of standard configuration for calibration use.
//
//    BrTpcTimeCalModule* t1timecalModule =
//      new BrTpcTimeCalModule("T1", "T1timerow.root");
//    t1timecalModule->SetMinFitHits(5);
//    t1timecalModule->SetSaveAscii(kTRUE);
//    t1timecalModule->SetCalibFile(Form("t1timecalib%d.txt", 
//  				     runOption->GetValue()));
//    t1timecalModule->SetTreeOn(kTRUE); 
//    t1timecalModule->AddFitRow( 3, 0); <-----|
//    t1timecalModule->AddFitRow( 4, 1); <---| |
//    t1timecalModule->AddFitRow( 7, 2); <-| | |
//    t1timecalModule->AddFitRow( 8, 3); <-------- Assumed good rows 
//    t1timecalModule->AddFitRow(11, 4); <-|  |    (row, index in array)
//    t1timecalModule->AddFitRow(12, 5); <----|
//    switcher->AddModule(t1timecalModule);
//

//____________________________________________________________________
//
// $Id: BrTpcTimeCalModule.cxx,v 1.3 2002/01/03 19:51:17 cholm Exp $
// $Author: cholm $
// $Date: 2002/01/03 19:51:17 $
// $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 BRAT_BrTableNames
#include "BrTableNames.h"
#endif
#ifndef BRAT_BrTpcTimeCalModule
#include "BrTpcTimeCalModule.h"
#endif
#ifndef ROOT_TF1
#include "TF1.h"
#endif
#ifndef ROOT_TProfile
#include "TProfile.h"
#endif
#ifndef ROOT_TNtuple
#include "TNtuple.h"
#endif
#ifndef ROOT_TString
#include "TString.h"
#endif
#ifndef BRAT_BrDataTable
#include "BrDataTable.h"
#endif
#ifndef BRAT_BrRunInfoManager
#include "BrRunInfoManager.h"
#endif
#ifndef BRAT_BrCalibration
#include "BrCalibration.h"
#endif
#ifndef BRAT_BrException
#include "BrException.h"
#endif

//____________________________________________________________________

ClassImp(BrTpcTimeCalModule);

//____________________________________________________________________
 BrTpcTimeCalModule::BrTpcTimeCalModule()
  : BrTpcCalModule(), kMaxRows(12)
{
  // Default constructor. DO NOT USE
  SetState(kSetup);

  fHitContainer = 0;
  fFitList = 0;

  SetMinFitHits();
  SetChi2Cut();
  SetDxCut();
}

//____________________________________________________________________
 BrTpcTimeCalModule::BrTpcTimeCalModule(const Char_t* name, 
						 const Char_t* title)
  : BrTpcCalModule(name, title), kMaxRows(12)
{
  // Named Constructor
  SetState(kSetup);

  fHitContainer = new TObjArray();
  fFitList = new TArrayI(kMaxRows);
  for(Int_t i = 0; i < kMaxRows; i++)  
    fFitList->AddAt(-1, i);
  
  SetMinFitHits();
  SetChi2Cut();
  SetDxCut();
}

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

  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; 
  
  fHistDir = gDirectory->mkdir(Form("%sRowTime",GetName())); 
  fHistDir->cd();
  
  hChi2    = new TH1F(Form("chi2%s", GetName()),
		      Form("%s:Chi2 for all hits", GetName()),
		      100, 0, 100);
  hChi2Acc = new TH1F(Form("chi2acc%s", GetName()),
		      Form("%s:Chi2 accepted", GetName()),
		      100, 0, 100);

  const Int_t nRows = fParamsTpc->GetNumberOfRows();

  hRowTime = new TH2F*[nRows];
  for(Int_t i = 0; i < nRows; i++)
    hRowTime[i] = 
      new TH2F(Form("hTpcDyVsY_row%d", i+1), 
	       Form("%s row time corrections. Row %d", 
		    GetName(), i+1),
	       100, 0, 160,
	       100, -5, 5);
  
  if(TreeOn()) {
    
    fTree = new TTree("T","Root tree");
    fTree->Branch("eventinfo", &fEventInfo, 
		  "x:y:dx:dy:chi2/F:row/I:hits/I:status/I");
  }

  // Histograms
  gDirectory = saveDir;
}

//____________________________________________________________________
 void BrTpcTimeCalModule::Init()
{
  // Job-level initialisation
  SetState(kInit);
  
  // base class initialization (register calibration parameters)
  BrTpcCalModule::Init();
  
  BrCalibration::EAccessMode mode = BrCalibration::kRead;
  
  const Int_t nRows = fParamsTpc->GetNumberOfRows();

  try{
    // check if we want to load adc gain cal from ascii file
    if (fLoadAscii) {
      mode = BrCalibration::kTransitional;
      fCalibration->Use(BrTpcCalibration::kTimeOffset0, mode, nRows);
      fCalibration->Use(BrTpcCalibration::kTimeOffset1, mode, nRows);
    } else {
      
      mode = BrCalibration::kWrite;
      
      fCalibration->Use(BrTpcCalibration::kTimeOffset0, mode, nRows);
      fCalibration->Use(BrTpcCalibration::kTimeOffset1, mode, nRows);
    }
  }
  catch(BrException* e) {
    cerr << "BrTpcTimeCalModule::Init : " << *e << endl;
    e->Execute();
  }
  
}


//____________________________________________________________________
 void BrTpcTimeCalModule::Begin()
{
  // Run-level initialisation
  SetState(kBegin);
  
  if (fLoadAscii) {
    ReadAscii();
    SetTimeOffsets();
    return;
  }
  
  if (!HistBooked())
    if (!fCommitAscii) {
      Abort("Begin", "Need histos for calibration");
      return;
    }
  
  // check if we got parameter revisions:
  if (!fCalibration->RevisionExists(BrTpcCalibration::kTimeOffset0) || 
      !fCalibration->RevisionExists(BrTpcCalibration::kTimeOffset1)) {
    Abort("Begin", 
	  "Some calibration revision are missing! BrDbUpdateModule ?");
    return;
  }
  
  // Set the corresponding parameters in BrDetectorParamsTPC to 0
  ZeroDetectorParameters();
}

//____________________________________________________________________
 void BrTpcTimeCalModule::ZeroDetectorParameters()
{
  // Set Row Time Corrections in BrDetectorParamsTPC
  const Int_t nRows = fParamsTpc->GetNumberOfRows();
  for(Int_t i = 0; i < nRows; i++) {
    const Int_t row = i+1;
    fParamsTpc->SetTimeCorrection(row, 0, 0);
  }

  fParamsTpc->ListParameters();
}

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

  if (fCommitAscii || fLoadAscii)
    return;
  
  // Get the data table with the tpc tracks
  BrDataTable *inputTable = 
    inNode->GetDataTable(Form("%s %s", BRTABLENAMES kTpcTrack, GetName()));
  
  if(!inputTable) {
    
    Warning("Event", "Table was not found : %s %s ", BRTABLENAMES kTpcTrack, GetName());
    return;
  }
  
  // loop over tpctracks
  const Int_t nTracks = inputTable->GetEntries();
  
  if(DebugLevel())
    cout << "Number of tracks : " << nTracks << endl;
  
  for(Int_t i = 0; i < nTracks; i++) {
    
    // for each track make a new special track where all the hit
    // positions have been converted to y positions and all errors set
    // to 1 Fit this new track with the middle rows only where the
    // time rows corrections should be smallest
    BrTpcTrackCandidate *tpcTrack = (BrTpcTrackCandidate*)inputTable->At(i);
    
    BrTpcTrackCandidate *newTrack = MakeNewTpcTrack(tpcTrack);
    
    // Fill the tree with the corrections for each row.
    if(newTrack) {
      
      FillHistograms(newTrack);
      
      delete newTrack;
      newTrack = 0;
    }

    // No matter what we have to delete the new hits
    fHitContainer->Delete();
  }
}

//____________________________________________________________________
 Bool_t BrTpcTimeCalModule::IsFitHit(BrTpcHit* hit)
{
  // If the row number is in fFitList return kTRUE
  // else return kFALSE
  const Int_t row = hit->GetRow();
  
  for(Int_t i = 0; i < kMaxRows; i++)
    if(row ==  fFitList->At(i)) {
      
      if(DebugLevel()>5)
	cout << "TRUE" << endl;

      return kTRUE;
    }
  
  if(DebugLevel()>5)
    cout << "FALSE" << endl;

  return kFALSE;
}

//____________________________________________________________________
BrTpcTrackCandidate* 
 BrTpcTimeCalModule::MakeNewTpcTrack(BrTpcTrackCandidate *track)
{
  // Loop over hits in track twice
  // First add all hits in fitrows and fit
  // Second add the rest of the hits 
  // If number of fit hits is lower than fMinFitHits return 0
  // else return new track
  if(DebugLevel()>2)
    cout << "In make new tpc track : " << endl;
  
  BrTpcTrackCandidate *newTrack = new BrTpcTrackCandidate();
  
  const Int_t nHits = track->GetNhit();
  
  for(Int_t i = 0; i < 2; i++) {
    
    for(Int_t j = 0; j < nHits; j++) {
      
      BrTpcHit* hit = track->GetTpcHitAt(j);
      
      if(DebugLevel()>5)
	cout << "I = " << i << "Checking for hit in row : " << hit->GetRow();
      
      if(i==0) {
	
	if(IsFitHit(hit))
	  AddTpcHit(newTrack, hit);
	
      } else if(i==1) {
	
	if(!IsFitHit(hit))
	  AddTpcHit(newTrack, hit);
      }
      
    }
    
    if(i==0) {
      
      if(newTrack->GetNhit() >= fMinFitHits) {
	
	newTrack->Fit();
      } else {

	delete newTrack;
	return 0;
      }
    }
  }
  
  return newTrack;
}

//____________________________________________________________________
 void BrTpcTimeCalModule::FillHistograms(BrTpcTrackCandidate* track)
{
  // Loop over tpc hits and fill histograms and tree, the last is optional
  if(DebugLevel()>1)
    cout << "In FillHistograms" << endl;
  
  const Int_t nHits = track->GetNhit();
  
  for(Int_t i = 0; i < nHits; i++) {
    
    const Float_t chi2 = track->GetQuality();
    
    // Get hitInfo
    BrTpcHit* hit = track->GetTpcHitAt(i);
    
    BrVector3D hitPos = hit->GetPos();
    
    //      // fill offsets
    BrVector3D fitPos;
    Float_t x, y;
    track->GetXYPositionAtZ(x, y, hitPos[2]);
    fitPos[0] = x;
    fitPos[1] = y;
    BrVector3D deviationVec = hitPos - fitPos;

    const Int_t row = hit->GetRow();

    if(TreeOn()) {
      
      if(DebugLevel()>5)
	cout << "Filling tree" << endl;

      fEventInfo.hits = nHits;
      fEventInfo.chi2 = chi2;
      fEventInfo.x = hitPos[0];
      fEventInfo.y = hitPos[1];
      fEventInfo.row = row;
      fEventInfo.status = hit->GetStatus();
      fEventInfo.dx = deviationVec[0];
      fEventInfo.dy = deviationVec[1];
      
      fTree->Fill();
    }
    
    hChi2->Fill(chi2);

    if(chi2 < fChi2Cut && TMath::Abs(deviationVec[0]) < fDxCut) {
      if(DebugLevel()>5)
	cout << "Filling histograms " << endl;
      hRowTime[row-1]->Fill(hitPos[1], deviationVec[1]);
      hChi2Acc->Fill(chi2);
    }
  }
}

//____________________________________________________________________
 void BrTpcTimeCalModule::AddTpcHit(BrTpcTrackCandidate* track, BrTpcHit* hit)
{
  // Add hit to track after converting to Time and Pad coordinates
  // Also store it in a container so we can delete it after use
  // 
  BrTpcHit* newHit = new BrTpcHit(*hit);
  
  const Float_t x = newHit->GetPos()[0];
  const Float_t y = newHit->GetPos()[1];
  newHit->SetX(fParamsTpc->IntrinsicToPad(x));
  newHit->SetY(fParamsTpc->IntrinsicToTime(newHit->GetRow(), y));
  
  const Double_t dx = TMath::Abs(fParamsTpc->IntrinsicToPad(0.05)-
				 fParamsTpc->IntrinsicToPad(0));
  const Double_t dy = TMath::Abs(fParamsTpc->IntrinsicToTime(0.05)-
				 fParamsTpc->IntrinsicToTime(0));
  newHit->SetXError(dx);
  newHit->SetYError(dy);

  track->AddHit(newHit);
  fHitContainer->Add(newHit);
}

//____________________________________________________________________
 void BrTpcTimeCalModule::Finish()
{
  // Job-level finalisation
  // This is very the fit to the histograms are done and the output
  // and commit is done
  SetState(kFinish);
  
  // if load ascii mode
  if (fLoadAscii) 
    return;
  
  // if commit mode
  if (fCommitAscii) {
    ReadAscii();
    return;
  }

  if (fSaveAscii) {

    cout << "In Finish" << endl;

    // ------------------------------------------------
    // fit 2d histograms with first degree polunomial
    // ------------------------------------------------
    
    TDirectory* saveDir = gDirectory; 
    
    fHistDir->cd();
    
    const Int_t nRows = fParamsTpc->GetNumberOfRows();

    for(Int_t j = 0; j < nRows; j++) {
      
      if(hRowTime[j]->GetEntries() == 0) {

	fCalibration->SetTimeOffset0(j+1, 0);
	fCalibration->SetTimeOffset1(j+1, 0);
	delete hRowTime[j];
	hRowTime[j] = 0;
	continue;
      }

      // Find the range in which to fit
      TH1D *tempHist = hRowTime[j]->ProjectionX();
      
      Double_t range[2];
      FindRanges(tempHist, range);
      delete tempHist;
    
      // Make a y profile of the 2d histogram
      // Remove the bins with few counts 
      // Fit with a first degree polunomial
      // Store results
    
      TProfile *prof = ProjectUsingGaus(hRowTime[j], 2);
    
      RemoveLowError(prof);
    
      prof->Fit("pol1","Q0","", range[0], range[1]);

      Char_t profName[64]; 
      sprintf(profName, "Profile:%s", hRowTime[j]->GetName());
      prof->SetName(profName);

      TF1 *func = prof->GetFunction("pol1");
      const Float_t timeBucket = fParamsTpc->GetTimeBucket();
      Double_t par0 = timeBucket * func->GetParameter(0);
      Double_t par1 = timeBucket * func->GetParameter(1);

      fCalibration->SetTimeOffset0(j+1, par0);
      fCalibration->SetTimeOffset1(j+1, par1);
    }

    // restore directory
    gDirectory = saveDir;

    // save calibration to ascii file
    SaveAscii();
  }
}

//_____________________________________________________________________
 void BrTpcTimeCalModule::RemoveLowError(TProfile *hist )
{
  // Remove bins in Profile histograms where there is less than 10
  // counts and the gaussian approximation of errors (sqrt(n)) is bad
  Int_t nBinX = hist->GetNbinsX();
  
  for(Int_t i = 1; i <= nBinX; i++) {
    if(hist->GetBinEntries(i) < 10){
        
        hist->SetBinEntries(i, 0);
        hist->SetBinContent(i, 0.0);
        hist->SetBinError(i, 0.0);
        
    } // end entries check
      
  } // end loop over x bins

}

//_____________________________________________________________________
 void BrTpcTimeCalModule::FindRanges(TH1D *hist, Double_t *range)
{
  // Find the range where we want to fit
  // low = first bin where n > n_max/3
  // high = last bin where n > n_max/3
  const Int_t nBinX = hist->GetNbinsX();
  const Double_t max = hist->GetMaximum();
  
  // Find min
  for(Int_t i = 1; i <= nBinX; i++) {
    if(hist->GetBinContent(i) > max/3) {
      range[0] = hist->GetBinLowEdge(i+1);
      break;
    }
  } 

  // Find max
  for(Int_t i = nBinX; i >= 1; i--) {
    if(hist->GetBinContent(i) > max/3) {
      range[1] = hist->GetBinLowEdge(i);
      break;
    }
  } 

}

//_____________________________________________________________________
 TProfile* BrTpcTimeCalModule::ProjectUsingGaus(TH2F *hist, Int_t nBins)
{
  // Project the 2d histogram in small x bins (size nBins)
  // The distribution is fitted with a gauss
  Int_t nBinX = hist->GetNbinsX();
  Int_t ratio = nBinX/nBins;
  TProfile *result = new TProfile("result", "Result", ratio, 0, 160,-50, 50);
  for(Int_t i = 1; i <= ratio; i++) {
    
    Int_t firstBin = nBins*(i-1) + 1;
    Int_t lastBin  = nBins*i;
    
    TH1D* projY = hist->ProjectionY("dummy", firstBin, lastBin);
    projY->Fit("gaus","Q0+");
    TF1 *func = projY->GetFunction("gaus");
    Float_t mean = func->GetParameter(1);
    //    Float_t sigma = func->GetParameter(2);
    Float_t sigma = func->GetParError(1);
    
    result->SetBinEntries(i, projY->GetEntries());
    result->SetBinContent(i, mean*projY->GetEntries());
    result->SetBinError(i, sigma*projY->GetEntries());
    delete projY;
  } 

  return result;
}

//____________________________________________________________________
 void BrTpcTimeCalModule::SaveAscii() 
{
  // save pedestal to ascii file
  
  cout << "In SaveAscii" << endl;

  BrRunInfoManager* runMan = BrRunInfoManager::Instance();
  const BrRunInfo*     run = runMan->GetCurrentRun();
  if (run->GetRunNo() == -1) {
    Stop("SaveAscii", "RunInfo has run number = -1");
    return;
  }
  
  cout << "Opening file" << endl;

  BrTpcCalModule::SaveAscii();
  
  ofstream file(fCalibFileName.Data(), ios::out);
  
  file << "****************************************** " << endl;
  file << "*  Calibration for Tpc detector " << GetName() << endl;
  file << "*  Time calibration         " << endl;
  file << "*  Used events from run " << run->GetRunNo() << endl;
  file << "****************************************** " <<endl;
  file << "*" << endl;    
  file << "* slat  |     top    |    bot     " << endl;
  file << "* --------------------------------" << endl << endl;
  
  const Int_t nRows = fParamsTpc->GetNumberOfRows();
  
  for (Int_t i = 0; i < nRows; i++) {
    const Int_t row = i + 1;
    file << setw(4)  << row 
	 << setw(15) << fCalibration->GetTimeOffset0(row) 
	 << setw(15) << fCalibration->GetTimeOffset1(row) << endl;
  }
  
  file << "* ------------------------------------" << endl << endl;
  
}


//____________________________________________________________________
 void BrTpcTimeCalModule::ReadAscii() 
{
  // read calibration from file created by this module

  BrTpcCalModule::ReadAscii();
  
  ifstream file(fCalibFileName.Data(), ios::in);
  
  if (!file) {
    Stop("ReadAscii", "File %s was not found", fCalibFileName.Data());
    return;
  }
  
  Char_t comment[256];
  
  file.getline(comment, 256);
  while(comment[0] == '*') {
    file.getline(comment, 256);
    if (DebugLevel() > 5)
      cout << comment << endl;
  } 

  Float_t timeCorrection0, timeCorrection1;
  Int_t   row;

  const Int_t nRows = fParamsTpc->GetNumberOfRows();
  
  for (Int_t i = 0; i < nRows; i++) {
    file >> row >> timeCorrection0 >> timeCorrection1;
    if (DebugLevel() > 5) 
      cout << setw(5)  << row 
	   << setw(12) << timeCorrection0 
	   << setw(12) << timeCorrection1 << endl;
    
    fCalibration->SetTimeOffset0(row, timeCorrection0);
    fCalibration->SetTimeOffset1(row, timeCorrection1);
  }
  
  fCalibration->SetComment(BrTpcCalibration::kTimeOffset0,
			   "Generated by BrTpcTimeCalModule: "
			   "pol1 fit to deviations");
  fCalibration->SetComment(BrTpcCalibration::kTimeOffset1,
			   "Generated by BrTpcTimeCalModule: "
			   "pol1 fit to deviations");
}

//____________________________________________________________________
 void BrTpcTimeCalModule::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: cholm $" << endl  
         << "    $Date: 2002/01/03 19:51:17 $"   << endl 
         << "    $Revision: 1.3 $ " << endl  
         << endl 
         << "-------------------------------------------------" << endl;
}

//____________________________________________________________________
//
// $Log: BrTpcTimeCalModule.cxx,v $
// Revision 1.3  2002/01/03 19:51:17  cholm
// Prepared to use BrTableNames class (or perhaps BrDetectorList) for table names
//
// Revision 1.2  2001/11/02 13:38:54  pchristi
// Added new module for calibrating drift velocities using the fibres behind
// T2 and in front of T1.
// Updated old modules with small changes.
//
// Revision 1.1  2001/10/12 11:08:50  pchristi
// Added new directory tpc. Added the first calibration modules. They
// have all been tested and found to work. The algorithms might not be
// optimal, but are at least fully automatic and to some extent
// documented.  CVS:
// ----------------------------------------------------------------------
// BrTpcCalModule.cxx BrTpcCalModule.h CVS: BrTpcHitWidthCalModule.cxx
// BrTpcHitWidthCalModule.h CVS: BrTpcPadStatusCalModule.cxx
// BrTpcPadStatusCalModule.h CVS: BrTpcTimeCalModule.cxx
// BrTpcTimeCalModule.h Include.h CVS: LinkDef.h Makefile.am CVS:
// ----------------------------------------------------------------------
//
// Revision 1.3  2001/10/10 15:49:55  pchristi
// Many updates and a new module
//
// Revision 1.2  2001/10/08 10:30:28  pchristi
// Classes compiles now
//
// Revision 1.1  2001/10/08 08:08:16  pchristi
// Initial revision.
//
//
//

This page automatically generated by script docBrat by Christian Holm

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

Validate HTML
Validate CSS