|
//____________________________________________________________________ // // Module to do the find the centrality cuts for TMA and SMA. // // The cuts are defined in terms of the total energy deposited rather // than the total multiplicity, so the the cuts (calibrations) are // largly indenpendent of the various corrections one may apply to // the multiplicty calculations. // // Cutting in the energy rather than the multiplicity ofcourse gives // something different, but I believe (and also have seen) that the // effect is negiable. Also, cutting in energy is less likely to // introduce system errors or bias. // // This means that we define the centrality in terms of energy // deposited in the TMA/SMA and _not_ in terms of multiplicity. // // The module expects as input minimum bias events. For the summer // 2000 run this is defined as // // Inclusive trigger 4 // At least one hit above thresshold in TMA/SMA // The existence of a primary vertex as determind by the ZDCs // // An additional vertex cut of +/- 50cm is highly recommended. Note // that only the primary vertex position as found from the ZDC data // is used, to avoid introducing a further bias. // //____________________________________________________________________ // // $Id: BrMultCentCalModule.cxx,v 1.2 2001/11/21 19:13:30 cholm Exp $ // $Author: cholm $ // $Date: 2001/11/21 19:13:30 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #ifndef BRAT_BrMultCentCalModule #include "BrMultCentCalModule.h" #endif #ifndef BRAT_BrDetectorList #include "BrDetectorList.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 ROOT_TGraph #include "TGraph.h" #endif #ifndef ROOT_TF1 #include "TF1.h" #endif #ifndef WIN32 #ifndef __IOSTREAM__ #include <iostream> #endif #ifndef __FSTREAM__ #include <fstream> #endif #ifndef __IOMANIP__ #include <iomanip> #endif #ifndef __CLIMITS__ #include <climits> #endif #ifndef __CFLOAT__ #include <cfloat> #endif #else #include <iostream.h> #include <fstream.h> #include <iomanip.h> #include <limits.h> #include <float.h> #endif //____________________________________________________________________ ClassImp(BrMultCentCalModule); //____________________________________________________________________ BrMultCentCalModule::BrMultCentCalModule() { // Default constructor. DO NOT USE SetState(kSetup); SetMax(); SetSingleNBins(); SetVertexMax(); SetVertexNBins(); SetFunctionOrder(); fNCuts = 0; fCuts = 0; fSCuts = 0; } //____________________________________________________________________ BrMultCentCalModule::BrMultCentCalModule(const Char_t* name, const Char_t* title) : BrModule(name, title) { // Named Constructor SetState(kSetup); SetMax(); SetSingleNBins(); SetVertexMax(); SetVertexNBins(); SetFunctionOrder(); fNCuts = 0; fSCuts = 20; fCuts = new Float_t[fSCuts]; SetOutputName(); fVtxVsEHisto = 0; fCutHistos = 0; } //____________________________________________________________________ BrMultCentCalModule::~BrMultCentCalModule() { // Destructor. if (!HistBooked()) return; } //____________________________________________________________________ void BrMultCentCalModule::AddCut(Float_t cut) { // Add a centrality cut to the calibration. Cut must be between 0 // (inclusive) and 100 (exclusive). if (GetState() != kSetup) { Failure("AddCut", "must be send before Init"); return; } // Warn if the cut is outside of range if (cut < 0 || cut >= 100) { Warning("AddCut", "'%f' invalid cut value, cut must be in [0,100)", cut); return; } // Expand the array if needed. if (fNCuts == fSCuts) { fSCuts = Ssiz_t(1.5 * fSCuts); Float_t* cuts = new Float_t[fSCuts]; for (Int_t i = 0; i < fNCuts; i++) cuts[i] = fCuts[i]; delete [] fCuts; fCuts = cuts; } fCuts[fNCuts++] = cut; } //____________________________________________________________________ void BrMultCentCalModule::SetFunctionOrder(Int_t order) { // Set the order of the polynomial for the (vertex,energy), // fit. Default is 3. if (GetState() != kSetup) { Failure("SetFunctionOrder", "must be send before Init"); return; } fFunctionOrder = order; } //____________________________________________________________________ void BrMultCentCalModule::SetOutputName(const Char_t* name) { // Set the file name of the file to output temporary calibration // parameters to if (GetState() != kSetup) { Failure("SetOutputName", "must be send before Init"); return; } fOutputName = name; } //____________________________________________________________________ void BrMultCentCalModule::SetMax(Float_t max) { // Set the maximum value to put in the internal 2D histogram for the // sum total calibrated tile response. Default is 4 if (GetState() != kSetup) { Failure("SetMax", "must be send before Init"); return; } fMax = max; } //____________________________________________________________________ void BrMultCentCalModule::SetSingleNBins(Int_t bins) { // Set the number of single bins. Default is 100. if (GetState() != kSetup) { Failure("SetSingleNBins", "must be send before Init"); return; } fSingleNBins = bins; } //____________________________________________________________________ void BrMultCentCalModule::SetVertexMax(Float_t max) { // Set the maximum Z-coordinate of the primary verticies to use in // fits. Default is 50cm. if (GetState() != kSetup) { Failure("SetVertexMax", "must be send before Init"); return; } fVertexMax = max; } //____________________________________________________________________ void BrMultCentCalModule::SetVertexNBins(Int_t bins) { // Set the number of vertex bins to fit over. // Default is 8. if (GetState() != kSetup) { Failure("SetVertexNBins", "must be send before Init"); return; } fVertexNBins = bins; } //____________________________________________________________________ void BrMultCentCalModule::DefineHistograms() { // Define histograms. They are: // <fill in here> if (GetState() != kInit) { Stop("DefineHistograms", "Must be called after Init"); return; } TString dirName(GetName()); dirName += "CentCalib"; TDirectory* saveDir = gDirectory; fHistDir = gDirectory->mkdir(dirName.Data()); if (!fHistDir) { Warning("DefineHistograms","could not create histogram subdirectory"); return; } fHistDir->cd(); // Make histograms here // Histograms for energy based cuts. fVtxVsEHisto = new TH2F("vtxVsE", "#frac{d^{2}N}{dV_{z}dE}", fVertexNBins, -fVertexMax, fVertexMax, fSingleNBins, 0, fMax); fVtxVsEHisto->SetXTitle("V_{z} [cm]"); fVtxVsEHisto->SetYTitle("E [?]"); fVtxVsEHisto->SetZTitle("#frac{d^{2}N}{dV_{z}dE}"); fVtxHisto = new TH1D("vtx", "Vertex distribution", fVertexNBins, -fVertexMax, fVertexMax); // Histograms for energy based cuts. fCutHistos = new TH1D * [fNCuts]; for (Int_t i = 0; i < fNCuts; i++) fCutHistos[i] = 0; for (Int_t i = 0; i < fNCuts; i++) { // Histograms for energy based cuts. fCutHistos[i] = new TH1D(Form("cut%02d", Int_t(fCuts[i])), Form("%02d%% Centrality Cut as function of Vz", Int_t(fCuts[i])), fVertexNBins, -fVertexMax, fVertexMax); fCutHistos[i]->SetXTitle("V_{z} [cm]"); fCutHistos[i]->SetYTitle("E_{cut}"); } gDirectory = saveDir; } //____________________________________________________________________ void BrMultCentCalModule::Init() { // Initialise this module. // // If no centrality cuts have been defined by user, then the // following cuts are defined: // // 5%, 10%, 15%, 20%, 30%, 40%, 50% // if (fNCuts == 0) { AddCut(5); // 5% centrality cut AddCut(10); // 10% centrality cut AddCut(15); // 15% centrality cut AddCut(20); // 20% centrality cut AddCut(30); // 30% centrality cut AddCut(40); // 40% centrality cut AddCut(50); // 50% centrality cut } // Must be here SetState(kInit); } //____________________________________________________________________ void BrMultCentCalModule::Begin() { // Run-level initialisation SetState(kBegin); if (!HistBooked()) Failure("Begin", "must book histograms for calibration to work"); } //____________________________________________________________________ void BrMultCentCalModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // Per event method SetState(kEvent); // Fill the internal histogram (vertex Z, tile response) SetState(kEvent); // Only go on if histograms are on if (!HistOn()) return; TString tableName; tableName = "Rdo"; tableName += GetName(); BrMultRdo* rdo = (BrMultRdo*)inNode->GetObject(tableName.Data()); if(!rdo) { Warning("Event", "couldn't get RDO for %s '%s'", GetName(), tableName.Data()); return; } // Try to get a vertex. Float_t vertexZ = FindVertex(inNode); if (vertexZ == FLT_MAX) vertexZ = FindVertex(outNode); if (vertexZ == FLT_MAX) { Stop("Event", "Cannot find z coordinate of primary vertex in node"); return; } fVtxHisto->Fill(vertexZ); // We require at least on hit over pedestal // if (rdo->GetArrayMultiplicity() < 1) if (rdo->GetArrayHits() < 1) { if (DebugLevel() > 0) cout << "<BrMultCentCalModule::Event>: No hits" << endl; return; } Float_t vtxPitch = 2 * fVertexMax / fVertexNBins; if (DebugLevel() > 0) cout << "<BrMultCentCalModule::Event>: Filling scatter plot " << vertexZ << ", " << rdo->GetArrayEnergy() << endl; fVtxVsEHisto->Fill(vertexZ, rdo->GetArrayEnergy()); } //__________________________________________________________________ Float_t BrMultCentCalModule::FindVertex(BrEventNode* node) { // Try to get a vertex from the input node. Only ZDC vertex is // considered to not introduce further bias // Look for a ZDC vertex TString tableName = "Rdo"; tableName += BrDetectorList::GetDetectorName(kBrZDC); BrZdcRdo* rdoZdc = (BrZdcRdo*)node->GetObject(tableName.Data()); if (rdoZdc && rdoZdc->GetZ() < 9999) return rdoZdc->GetZ(); return FLT_MAX; } //____________________________________________________________________ void BrMultCentCalModule::Finish() { // Do the calibration. That is, do the linear (polynomial) fit to // the values of the cuts requested. The parameters is written to a // ASCII file on disk, and histograms are stored in reqular // histogram file. SetState(kFinish); if (DebugLevel() > 4) cout << "<BrMultCentCalModule::Finish()> for " << GetName() << endl; // We need the histograms if (!HistBooked()) { Failure("Finish", "must book histograms for calibration to work"); return; } TDirectory* savDir = gDirectory; fHistDir->cd(); // Try to open the output file ofstream outputFile(fOutputName.Data(), ios::out|ios::trunc); if(!outputFile) { Failure("Finish", "couldn't open file '%s'", fOutputName.Data()); return; } // Needed for Microshit Visual F**k++ and KCC Int_t i = 0, j = 0, k = 0; // Distance between vertex cuts. Float_t vtxPitch = 2 * fVertexMax / fVertexNBins; // Loop over vertex bins for (i = 0; i < fVertexNBins; i++) { // Integral of slice at constant vertex Z (one bin) Stat_t full = fVtxVsEHisto->Integral(i, i+1, 0, fSingleNBins); Short_t cutBin[fNCuts]; for (j = 0; j < fNCuts; j++) cutBin[j] = -1; Float_t oldratio = 0; // Loop over upper edge for intergal for (j = fSingleNBins - 1; j > 0; j--) { Stat_t integral = fVtxVsEHisto->Integral(i, i+1, 0, j); Float_t ratio = (1 - Float_t(integral) / full) * 100; // Loop over cuts for (k = 0; k < fNCuts; k++) { // We've found a bin that is higher fraction of all events // then cut, so set the cut bin to be the previous bin. But // only do it first time we are are above cut. if (fCuts[k] < ratio && cutBin[k] == -1) { Double_t binx = (i + .5) * vtxPitch - fVertexMax; Double_t biny = fVtxVsEHisto->GetYaxis()->GetBinCenter(j); fCutHistos[k]->Fill(binx, biny); cutBin[k] = j; } } // Done with cuts at this bin } // Done with cuts at this vertex bin // delete [] cutBin; } // Done with all vertex bins // Write a header outputFile << "# N single bins: " << fSingleNBins << endl; outputFile << "# N vertex bins: " << fVertexNBins << endl; outputFile << "Order = " << fFunctionOrder << "; " << endl; outputFile << "NCuts = " << fNCuts << "; " << endl; // Header for the fits outputFile << "# Cut functions" << endl << "# Parameters" << endl << "# Centrality " << flush; for (i = 0; i <= fFunctionOrder; i++) outputFile << setw(12) << i << " " << flush; outputFile << endl; // Now it's time to fit // TF1::InitStandardFunctions(); for (i = 0; i < fNCuts; i++) { // Make the fit, and store function if (DebugLevel() > 9) cout << GetName() << " Cent Calib: Fitting " << endl; TString e(Form("pol%d",fFunctionOrder)); // Char_t e[6]; // sprintf(e, ; fCutHistos[i]->Fit(e.Data(),"Q0+"); TF1* f1 = fCutHistos[i]->GetFunction(e); if (DebugLevel() > 9) { cout << GetName() << " Result of fit: " << flush; f1->Print(); } // Write the centrality first outputFile << setw(12) << fCuts[i] << " " << flush; // Write out parameters of fit. for (j = 0; j <= fFunctionOrder; j++) outputFile << setw(12) << f1->GetParameter(j) << " " << flush; outputFile << endl; } // Fitting loop if (DebugLevel() > 9) cout << GetName() << " Cent Calib: all done" << endl; if (DebugLevel() > 5) gDirectory->ls(); savDir->cd(); outputFile.close(); } //____________________________________________________________________ void BrMultCentCalModule::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: Christian Holm Christensen" << endl << " Last Modifications: " << endl << " $Author: cholm $" << endl << " $Date: 2001/11/21 19:13:30 $" << endl << " $Revision: 1.2 $ " << endl << endl << " Centrality cuts to determine:" << flush; for (Int_t i = 0; i < fNCuts; i++) cout << setw(5) << fCuts[i] << flush; cout << endl << " Polynomial order: " << fFunctionOrder << endl << " Vertex bins: " << fVertexNBins << endl << " Max absolute vertex Z: " << fVertexMax << endl << " Max signal: " << fMax << endl << " Signal bins: " << fSingleNBins << endl << endl << " Output file name: " << fOutputName << endl << "-------------------------------------------------" << endl; } //____________________________________________________________________ // // $Log: BrMultCentCalModule.cxx,v $ // Revision 1.2 2001/11/21 19:13:30 cholm // Some fixes, and removed old unused and redundant code. // // Revision 1.1 2001/09/04 16:53:00 cholm // New centrality calibration module BrMultCentCalModule, based on cuts // in energy deposited in the single detectors SMA and TMA, as outline in // BAN 34. // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|