//____________________________________________________________________
//
// Determine the centrality of the event from the full multiplicity
// array Reduced data objcts (Rdo). Input event must contain BrMultRdo
// and BrMultRdo objects. The primary vertex Z position need to be
// known, so either BrZdcRdo, BrBbRdo, or BrVertex from TPM1 should be
// present in the inout node. 
//
 
//
// $Id: BrMultCentModule.cxx,v 1.15 2002/03/23 16:09:12 sanders Exp $ 
// $Author: sanders $  
// $Date: 2002/03/23 16:09:12 $
// $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov>
//
#ifndef BRAT_BrMultCentModule
#include "BrMultCentModule.h"
#endif
#ifndef BRAT_BrTableNames
#include "BrTableNames.h"
#endif
#ifndef BRAT_BrDetectorList
#include "BrDetectorList.h"
#endif
#ifndef BRAT_BrMultCent
#include "BrMultCent.h"
#endif
#ifndef BRAT_BrMultRdo
#include "BrMultRdo.h"
#endif
#ifndef BRAT_BrVertex
#include "BrVertex.h"
#endif
#ifndef BRAT_BrBbRdo
#include "BrBbRdo.h"
#endif
#ifndef BRAT_BrZdcRdo
#include "BrZdcRdo.h"
#endif
#ifndef BRAT_BrParameterDbManager
#include "BrParameterDbManager.h"
#endif
#ifndef BRAT_BrCalibrationManager
#include "BrCalibrationManager.h"
#endif
#ifndef BRAT_BrMultCentCalibration
#include "BrMultCentCalibration.h"
#endif
#ifndef ROOT_TString
#include "TString.h"
#endif
#ifndef ROOT_TMath
#include "TMath.h"
#endif
#ifndef ROOT_TDirectory
#include "TDirectory.h"
#endif
#ifdef WIN32
#include <limits.h>
#include <float.h>
#include <iostream.h>
#include <iomanip.h>
#else 
#include <climits>
#include <cfloat>
#include <iostream>
#include <iomanip>
#endif
#define BR_CENT_MAX_VERTEX 35.

//____________________________________________________________________
ClassImp(BrMultCentModule);

//____________________________________________________________________
 BrMultCentModule::BrMultCentModule()
{
  // Default constructor
  SetSiWeight();
}

//____________________________________________________________________
 BrMultCentModule::BrMultCentModule(const Char_t* name, 
				   const Char_t* title)
  : BrModule(name, title), BrMultUtil(name)
{
  // Default constructor. 
  SetSiWeight();
}

//____________________________________________________________________
 void BrMultCentModule::DefineHistograms() 
{
  // Define the histograms for this module 
  // Currently, only a high vs low cut histogram is made. 

  // remember current directory
  TDirectory* savdir = gDirectory;   
  TDirectory* hdir   = savdir->mkdir("multcent");
  
  if (!hdir) {
    Warning("DefineHistograms","could not create histogram subdirectory");
    return;
  }
  
  hdir->cd();
  
  fCentHist = new TH1D("multcent", "Centrality", 102, 0, 101); 
  fMultHist = new TH1D("multAverage",
		       "Average Mult of Si and Tile",500,0,3000);
  fMultSiHist = new TH1D("multSi","Mult from Si",500,0,3000);
  fMultTileHist = new TH1D("multTile","Mult from Tile",500,0,3000);
  

  gDirectory = savdir;
  
}

//____________________________________________________________________
 void BrMultCentModule::Init() 
{
  // Initialise the object. Get the parameters for TMA or SMA, and
  // get a handle on the TMA or SMA calibrations.
  SetState(kInit);
  
#ifndef BR_MULT_CAL_TMP
  // Calibrations 
  BrCalibrationManager *manager = 
    BrCalibrationManager::Instance();
  fCalibration = static_cast<BrMultCentCalibration*>
    (manager->Register("BrMultCentCalibration", 
		       "Centrality"));
#else 
  fCalibration     = BrMultCentTmpCalibration::Instance();
  const BrRunInfo* runInfo = BrRunInfoManager::Instance()->GetCurrentRun();
  fCalibration->SetRunNumber( runInfo->GetRunNo() ); 
#endif 

  if (!fCalibration) {
    Failure("Init", "didn't get Centrality calibration from manager");  
    return;
  }

  fCalibration->Use("centralityFuncs");
}

//____________________________________________________________________
 void BrMultCentModule::Event(BrEventNode* input, BrEventNode* output) 
{
  // 
  Double_t multSi, multTiles, mult, cent;
  Double_t multDiff;
  
  SetState(kEvent);

  // Get the tile RDO 
  TString tableName;
  tableName = "Rdo";
  tableName += BrDetectorList::GetDetectorName(kBrTile);
  BrMultRdo* rdoTile = (BrMultRdo*)input->GetObject(tableName.Data());
  if(!rdoTile) {
    Warning("Event", "couldn't get Tile rdo's '%s'", 
	    tableName.Data());
    return;
  }
  // Get the Silicon RDO 
  tableName = "Rdo";
  tableName += BrDetectorList::GetDetectorName(kBrSi);
  BrMultRdo* rdoSi = (BrMultRdo*)input->GetObject(tableName.Data());
  if(!rdoSi) {
    Warning("Event", "couldn't get Si rdo's '%s'", 
	    tableName.Data());
    return;
  }

  // Make the output object. 
  tableName = "CentMult";
  BrMultCent* centobj = new BrMultCent(tableName.Data(),
			       "TMA + SMA centrality");
  output->AddObject(centobj);
  
  // Try to find the vertex. 
  Double_t vtx = FindVertex(input, output);

  // Get the multiplicites 
  multSi    = rdoSi->GetArrayMultiplicity();
  multTiles = rdoTile->GetArrayMultiplicity();
  //cout << "multSi = " << multSi << ", multTiles = " << multTiles << endl;
  mult      = -1000.;
  cent      = -1000.;

  const BrRunInfo* runInfo = BrRunInfoManager::Instance()->GetCurrentRun();
  if(runInfo && runInfo->GetRunNo() >= 3700) SetSiWeight(1.66); 
  
  // Check that we're inside the valid range of the calibrations. 
  if(TMath::Abs(vtx) < BR_CENT_MAX_VERTEX) {
    // Check that the multiplicites are ok 
    if (multSi > 0 && multTiles > 0)  
      // Calculate the weigted multiplicity 
      mult = (fSiWeight * multSi + multTiles) / 2.;
    // Calculate the weighted difference 
    multDiff = fSiWeight * multSi - multTiles;
    Float_t multDiffScaled = multDiff/(2.0*mult+10);
    // If the multiplicity is very low, or the diffence between the
    // two mulitplicities is 70%, relative to the summed multiplicity,
    // then calculate the centrality
    if (mult > 0.1 && multTiles>=4 &&
	TMath::Abs(multDiffScaled) < 0.8) 
      {
	cent = MultToCent(mult);
	if (HistOn()) {
	  fCentHist->Fill(cent);
	  fMultHist->Fill(mult);
	}
      }	
    if (HistOn()) {
      fMultSiHist->Fill(multSi);
      fMultTileHist->Fill(multTiles);
    }
  }
  

  // Set the output data 
  centobj->SetMult(mult);
  centobj->SetCent(cent);

  if (DebugLevel() > 10) 
    centobj->Print();
}

//____________________________________________________________________
 Double_t BrMultCentModule::MultToCent(Double_t mult)
{
  // PRIVATE METHOD:
  // Convert the multiplicity to centrality 
  Float_t cent;
  Double_t * params;
  Float_t x;
  Int_t order;
  Int_t i;
  x = mult;
  Bool_t earlyRun = kTRUE;

  // If the multiplicity is very low, return 100 automatically 
  if(mult<1) 
    return 100;

  // Cut of at high multiplicities.  This probably needs to be a
  // calibration constant, as the maximum of the multiplicity
  // distributiuon changes with square root s_NN 
  Float_t fMaxMult = 1680;
  const BrRunInfo* runInfo = BrRunInfoManager::Instance()->GetCurrentRun();
  if(runInfo && runInfo->GetRunNo() >= 3700)
    {
      fMaxMult = 2200;
      earlyRun = kFALSE;
    }
  if(mult>fMaxMult) 
    return 0.;
  //For 130 GeV runs we do not have a middle range
  if(earlyRun)
    { 
      // Else, if the multiplicity is high, use the Low parameterisation. 
      if(mult>100) {
	order = fCalibration->GetCentParamLOrder();
	params = fCalibration->GetCentParamL();
	cent = params[0] + params[1]*x;
	for(i=2;i<=order;i++)
	  cent+=params[i]*TMath::Power(x,i);
      }       
      // Else, if the multiplicity is low, use the High parameterisation. 
      else {
	order = fCalibration->GetCentParamHOrder();
	params = fCalibration->GetCentParamH();
	cent = params[0] + params[1]*x;
	for(i=2;i<=order;i++)
	  cent+=params[i]*TMath::Power(x,i);
      }
    }
  else
    {
      //Else, if the multiplicity is high, use the Low parameterisation. 
      if(mult>1950) 
      	{
	  order = fCalibration->GetCentParamLOrder();
	  params = fCalibration->GetCentParamL();
	  cent = params[0]*TMath::Exp(params[1]*(x-1900));;
	} 
      else if (mult>100)
	{      
	  // Else, if the multiplicity is mid, use the Mid parameterisation. 
	  order = fCalibration->GetCentParamMOrder();
	  params = fCalibration->GetCentParamM();
	  cent = params[0] + params[1]*x;
	  for(i=2;i<=order;i++)
	    cent+=params[i]*TMath::Power(x,i);
	}
      else 
	{
	  // Else, if the multiplicity is low, use the high parameterisation. 
	  order = fCalibration->GetCentParamHOrder();
	  params = fCalibration->GetCentParamH();
	  cent = params[0] + params[1]*x;
	  for(i=2;i<=order;i++)
	    cent+=params[i]*TMath::Power(x,i);
	}
    }
      
  // Make sure the centrality isn't lower than 0, or higher than 100
  if(cent < 0) 
    return 0;
  if(cent > 100) 
    return 100;

  // Return the thing. 
  return cent;
}

//____________________________________________________________________
//
// $Log: BrMultCentModule.cxx,v $
// Revision 1.15  2002/03/23 16:09:12  sanders
// Modified mutliplicity to centrality calibration function.  The
// high end of the mutliplicity distribution is now fitted with an
// exponential, rather than a high order polynomial. The maximum
// vertex for which the centrality calculation is performed has
// also been changed to 35 cm, rather than the previous 45 cm.
// It should be noted that even allowing 35 cm is dangerous.  The
// actual calibrations only use vertices < 30 cm....which is already
// outside of the si array.
//
// Revision 1.14  2002/02/21 23:28:49  sanders
// Restored vertex range to +/- 45 cm.  Also, reduced resolution on
// multiplicity histograms.
//
// Revision 1.13  2002/02/19 21:13:50  sanders
// added run number check for calibration
//
// Revision 1.12  2002/02/06 23:41:57  sanders
// Modified dispersion and range of histograms to make them more
// useful in generating the multiplicity to centrality calibration.
//
// Revision 1.11  2002/02/05 23:21:53  sanders
// Changed default vertex range to +-30 cm; added minimum tile mult
// cut to average mult.  These changes bring calculations in line with
// current centrality calibration procedure.
//
// Revision 1.10  2002/01/25 18:09:33  bjornhs
// Added a few more histograms
//
// Revision 1.9  2001/12/07 14:29:39  sanders
// another minor correction to check if -1 is added into a sum mulitiplicity
//
// Revision 1.8  2001/11/19 04:11:09  sanders
// Revised centrality calibration for 2001 (200 GeV) data.
//
// Revision 1.7  2001/11/13 15:51:16  cholm
// REmoved initialisation of calibrations.
//
// Revision 1.6  2001/11/12 14:57:03  sanders
// Added check for run info to allow use with 130 or 200 GeV data
// calibrations.
//
// Revision 1.5  2001/10/08 11:28:27  cholm
// Changed to use new DB access classes
//
// Revision 1.4  2001/09/20 13:14:48  cholm
// Changed to use BrMultUtil rather than BrMultModule.
//
// Revision 1.3  2001/09/04 16:58:01  cholm
// Minor clean up (removed some calibrations that wasn't used at all -
// very inefficient and potentially dangerous)
//
// Revision 1.2  2001/06/22 17:39:45  cholm
// Changes t odo with data classes renaming.
//
// Revision 1.1.1.1  2001/06/21 14:55:06  hagel
// Initial revision of brat2
//
// Revision 1.7  2001/06/02 15:23:26  sanders
// Modified parameters for E-to-N conversions.  Also modified parameters
// for converting the average Si+Tile sum multiplicity to centrality.  The
// new parameters result from an analysis where GEANT simulations are used
// to correct for the excessive number of low multiplicity events arising from
// pedestal tails being counted as particle hits.  Events with tile multiplicities
// between 4 and 200 are compared to the corresponding GBRAHMS (Hijing input)
// simulations to estimate the number of expected events with tile multiplicity
// below 4.
//
// Revision 1.6  2001/05/31 01:45:56  cholm
// Second rewamp of this directory.  All RDO modules use the common
// BrMultRdoModule, and they both write BrMultRdo Objects.  Also introduced
// ABC for Br[Tile|Si][Parameters|Calibration] since they have a lot in common,
// and the code can be made more efficient this way.
//
// A possible further thing to do, is to make an ABC for the CentModules, and
// the corresponding calibration modules, since they are very VERY similar.
// However, the current module BrMultCentModule is in the way.  Need to talk
// to Steve about that.
//
// Revision 1.5  2001/05/23 12:26:40  cholm
// Some important bug fixes, and proper histogram names
//
// Revision 1.4  2001/05/02 22:37:15  sanders
// Fixed minor mistake (that Hiro introduced??)
//
// Revision 1.3  2001/05/01 19:41:01  hito
// Minor bug fix
//
// Revision 1.2  2001/04/19 15:49:11  cholm
// Changed table name to correspond to the rest
//
// Revision 1.1  2001/04/19 15:09:50  cholm
// Renamed BrCent... to BrMultCent... to reduce confusion
//
//

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

Validate HTML
Validate CSS