/////////////////////////////////////////////////////////////////////////////
//
// EEezAnalysis
//
// Author: Jason C. Webb <jwebb@iucf.indiana.edu>
//
// The EEezAnalysis class is the main workhorse of the analysis code.
// It is responsible for taking the raw data from the ezTree, making
// the connection to the database to obtain and correct for pedestals
// and gains, and packing the data into the EEezTower objects. To keep
// the EEezTower class as simple as possible, all of the geometry
// information carried by the towers is initialized from the EEezAnalysis
// class. This has the side effect of making the cluster finder fairly
// independent of the EEMC geometry.
//
// The EEezAnalysis class has privately-held arrays to hold all detector
// elements: <br>
// EEezTower m_Towers[720];
// EEezStrip m_Strips[12][2][288];
//
// Several inline member functions are available which provide direct
// access to the m_Towers array through pointers. Each element in
// m_Towers contains pointers to neighboring towers, and to the
// appropriate strips in the m_Strips array.
//
// One accesses the data stored in the m_Towers array by calling member
// functions of EEezAnalysis. These member functions return pointers into
// the m_Towers array. One should tread lightly, as any changes you make
// to a tower will be made for all running code which uses this instance
// of EEezAnalysis. (This is why I have not told you about all of the
// "set" methods on towers and strips!)
//
// NOTE: One may place this class into an StMaker and use it to analyse
// muDst's. An example of this may be found in
//
// StRoot/StEEmcPool/StMuEEmcClusterMaker
//
///////////////////////////////////////////////////////////////////////////////
#include "EEezAnalysis.h"
#include <assert.h>
#include <iostream>
#include <stdio.h>
#include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
#include "StEEmcUtil/EEmcSmdMap/EEmcSmdMap.h"
#include "StEEmcUtil/EEfeeRaw/EEfeeRawEvent.h"
#include "StEEmcUtil/EEfeeRaw/EEfeeDataBlock.h"
#include "StEEmcUtil/EEfeeRaw/EEname2Index.h"
#include "StEEmcUtil/EEfeeRaw/EEdims.h"
#include "StEEmcUtil/EEevent/EEeventDst.h"
#include "StEEmcUtil/EEevent/EEsectorDst.h"
#include "StEEmcUtil/EEevent/EEtwHitDst.h"
#include "StEEmcUtil/EEevent/EEsmdHitDst.h"
#include "StEEmcDbMaker/cstructs/eemcConstDB.hh"
#ifndef __ROOT__
#include "EEmcDb.h"
#else
//#include "StEEmcDbMaker/StEEmcDbMaker.h"
#endif
#include "StEEmcDbMaker/EEmcDbItem.h"
typedef EEmcDbItem EEmcDbItem_t;
#include <TClonesArray.h>
#include <algorithm> // for std::random_shuffle
#ifdef __ROOT__
#include "StMuDSTMaker/COMMON/StMuEmcCollection.h"
#endif
// Kludge until we have gains...
//$$$#define EEMC_IDEAL_GAINS
//$$$#define ESMD_IDEAL_GAINS
#define EPRE_IDEAL_GAINS
/*
ETA IDEAL SLOPE ERR RMS
01 18.938 -0.23695 0.002600 0.02634
02 20.769 -0.23955 0.003385 0.03587
03 22.650 -0.20206 0.002979 0.03668
04 24.575 -0.21338 0.003129 0.03345
05 26.539 -0.20125 0.002710 0.03930
06 28.514 -0.18615 0.002483 0.02601
07 30.493 -0.16956 0.002537 0.02345
08 32.473 -0.16139 0.002180 0.01877
09 34.438 -0.14518 0.002103 0.01586
10 36.387 -0.13228 0.002369 0.02065
11 38.289 -0.12463 0.002141 0.01569
12 40.146 -0.12653 0.001870 0.02284
*/
Float_t ideal_gains[] = { 18.938, 20.769, 22.650, 24.575,
26.539, 28.514, 30.493, 32.473,
34.438, 36.387, 38.289, 40.146 };
Float_t ideal_smd = 23000;
Float_t ideal_pre = 23000;
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
ClassImp(EEezAnalysis);
EEezAnalysis::EEezAnalysis()
{
// Constructor
m_EEmcDb = NULL; // Start w/ NULL db
m_EEmcSmdMap = EEmcSmdMap::instance();
setSeedEnergy();
//$$$ m_Random = new TRandom(); // DANGER
m_StatMask = 0x0000;
// Initialize the name-to-index map for easy lookup of towers by name
for ( Int_t i = 0; i < 720; i++ ) {
Int_t iphi = ( i / 12 );
Int_t ieta = ( i % 12 ) + 1;
Int_t isec = ( iphi / 5 ) + 1;
Int_t isub = ( iphi % 5 ) + (Int_t)'A';
TString tow = "";
if ( isec < 10 ) tow += "0";
tow += isec;
tow += "T";
tow += (Char_t)isub;
if ( ieta < 10 ) tow += "0";
tow += ieta;
//$$$ std::cout << "Initializing tower " << i << " is " << tow << std::endl;
m_IndexFromName[ std::string( tow.Data() ) ] = i;
}
#ifndef __ROOT__
TString name = "EEsectorStats";
for ( Int_t i = 0; i < 12; i++ ) {
TString nname = name; nname += (i+1);
m_SectorStats[i] = new EEezSectorStats( nname.Data() );
}
#endif
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::addMask ( unsigned bit )
{
// Eliminates tower/strip/pre/postshower detectors which have
// the specified bit set in the database. The following bits
// have been defined in StEEmcDbMaker/cstructs/eemcConstDB.hh
//
// #define EEMCSTAT_ONLPED 0x0001 // only pedestal
// #define EEMCSTAT_STKBT 0x0002 // sticky lower bits
// #define EEMCSTAT_HOTHT 0x0004 // hot for HT trigger
// #define EEMCSTAT_HOTJP 0x0008 // hot for JP trigger
// #define EEMCSTAT_OUTPI0 0x0010 // hot in pi0 analysis
// #define EEMCSTAT_HOTSTR 0x0020 // hot esmd strip
//
// Once a mask has been set, it cannot be unset in the same analysis.
//
// Note that any detector with a nonzero "fail" bit will
// always be eliminated.
m_StatMask |= bit;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
EEezTower *EEezAnalysis::getTower ( Int_t index )
{
// Given the index (0-719) of an EEmc Tower, return a pointer to
// the EEezTower within the internal array.
assert ( 0 <= index && index < 720 );
return &m_Towers[index];
}
EEezTower *EEezAnalysis::getTower ( Int_t sector, Int_t subsector, Int_t etabin )
{
// Given the sector (0-11), subsector (0-4) and etabin (0-11), return
// a pointer to the appropriate EEezTower within the internal array.
return getTower( 60*sector + 12*subsector + etabin );
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::forceTower( const Char_t *tow, Float_t etow, Float_t epre1, Float_t epre2, Float_t epost )
{
// A bad idea, will probably delete in near future.
Int_t ii = m_IndexFromName[ std::string(tow) ];
assert ( 0 <= ii && ii < 720 );
m_Towers[ii].setEnergy(etow,0);
m_Towers[ii].setEnergy(epre1,1);
m_Towers[ii].setEnergy(epre2,2);
m_Towers[ii].setEnergy(epost,3);
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::Clear( Option_t *opts )
{
// Clear out towers and strips... as an optimization, we
// only clear out those towers and strips which have been
// hit.
EEezTowerPtrVecIter_t towerIter = m_HitTowers.begin();
while ( towerIter != m_HitTowers.end() ) {
(*towerIter)->clear();
towerIter++;
}
EEezStripPtrVecIter_t stripIter;
stripIter = m_HitUStrips.begin();
while ( stripIter != m_HitUStrips.end() ) {
(*stripIter)->clear();
stripIter++;
}
stripIter = m_HitVStrips.begin();
while ( stripIter != m_HitVStrips.end() ) {
(*stripIter)->clear();
stripIter++;
}
m_HitTowers.clear();
m_SeedTowers.clear();
m_HitUStrips.clear();
m_HitVStrips.clear();
#ifndef __ROOT__
for ( Int_t s = 0; s < 12; s++ ) { m_SectorStats[s] -> Clear(); }
#endif
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
#ifndef __ROOT__
void EEezAnalysis::Fill( EEfeeRawEvent *eevent )
{
// Takes as input a EEfeeRawEvent (most likely a branch from
// an ezTree), and fills the m_Towers and m_Strips arrays
// with ADC and reconstructed energy response.
// --
//
// Unpacks the raw fee event and fills a 720 element array
// of EEezTowers.
//
assert(m_EEmcDb); // Must setDb(EEmcDb *db) first
assert(eevent);
//
// Loop over all data blocks (i.e. crates and MAPMT boxen)
//
for ( Int_t ic = 0; ic < eevent -> block -> GetEntries(); ic++ ) {
//
// Get the data block from the event
//
EEfeeDataBlock *b = (EEfeeDataBlock *)eevent -> block -> At(ic);
// Validate
if( !b->isValid() ) continue;
// Handle incorrect crate id set on crate 1
if ( b -> getCrateID() == 0x10 ) b-> setCrateID(1);
//
// Get the data from the data block
//
Int_t crateID = b -> getCrateID();
UShort_t *data = b -> getData();
Int_t nd = b -> getValidDataLen();
//
// Loop over all channels in this crate
//
for ( Int_t chan = 0; chan < nd; chan++ ) {
//
// Access the database for information on gains, pedestals,
// pixel name, etc...
//
const EEmcDbItem_t *x = m_EEmcDb -> getByCrate(crateID,chan);
// Null pointer from DB
if ( x==0 ) continue;
// Channel is marked as bad in DB
if ( x -> fail ) continue;
//
// Mask off any channel which the user has requested
// to be masked off.
//
if ( m_StatMask & x -> stat ) continue;
Float_t value=data[chan]; // raw ADC
//--
//
// Handle towers, pre- and postshower readouts here
//
if ( !x -> isSMD() ) {
// Check if tower is marked as "hot", and
// skip it if it is.
if ( m_HotTowers.find( std::string( x->name ) ) !=
m_HotTowers.end() ) {
continue;
}
//
// Determine which sector, subsector, etabin was
// struck, what the pedestal and threshold is
// etc...
//
Int_t isec = (Int_t) (x -> sec - 1);
Int_t isub = (Int_t) (x -> sub - (Int_t)'A');
Int_t ieta = (Int_t) (x -> eta - 1);
//Int_t iphi = (Int_t) (5*isec + isub);
Int_t itow = 60*isec + 12*isub + ieta;
Float_t threshold = x -> thr;
//
// If this tower's ADC response is larger than
// the current "high-tower", set a new high-tower
//
if ( m_HighTower -> getRawAdc() > value )
m_HighTower = &m_Towers[itow];
//
// ADC must exceed the threshold specified w/in
// the database. If not, punt.
//
if ( value < threshold ) continue;
//
// Select which detector readout we have. Towers,
// pre and postshower, in that order.
//
Int_t idet = -1;
if ( x->name[2] == 'T' ) idet = 0;
else if ( x->name[2] == 'P' ) idet = 1;
else if ( x->name[2] == 'Q' ) idet = 2;
else if ( x->name[2] == 'R' ) idet = 3;
else
assert(0); // Houston, we have a problem.
#if 0 // Possible debugging code here
//Float_t ped = x -> ped;
std::cout << x -> name << std::endl;
std::cout << "isec = " << isec << std::endl
<< "isub = " << isub << std::endl
<< "ieta = " << ieta << std::endl;
// << "iphi = " << iphi << std::endl;
std::cout << "idet = " << idet << std::endl;
std::cout << "adc = " << value << std::endl;
std::cout << std::flush;
#endif
//
// Correct for gain to obtain the energy... Gains
// below zero get excluded from further analysis
// for obvious reasons.
//
#ifndef EEMC_IDEAL_GAINS
if ( x -> gain <= 0. && idet == 0 ) continue; // skip tower if db says to anyway
Float_t gain = x -> gain;
Float_t ped = x -> ped;
#endif
#ifdef EEMC_IDEAL_GAINS
Float_t gain = ideal_gains[ieta];
if ( idet != 0 ) gain = x -> gain;
Float_t ped = x -> ped;//dumb, for MC
#endif
#ifdef EPRE_IDEAL_GAINS
if ( idet > 0 ) gain = 23000; // see StEEmcFastSimulator
#endif
// If we've navigated the above and come up
// with a nonzero gain, we can process the
// event
if ( gain <= 0 ) continue;
#if 0
std::cout << "ped = " << ped << std::endl;
std::cout << "thr = " << threshold << std::endl;
std::cout << std::flush;
#endif
Float_t energy = ( value - ped ) / gain;
#if 0
std::cout << "ener = " << energy << std::endl;
std::cout << std::flush;
#endif
// if ( TString(x->name).Contains("08U") ||
// TString(x->name).Contains("08V") ) x->print();
//
// Punt if value < 0.0... this would indicate a
// negative gain (error condition) and should
// therefore be ignored.
//
if ( energy < 0.0 ) continue;
//
// Set the energy-response for the tower (pre or
// postshower). .
//
m_Towers[itow].setAdc( value-ped, idet );
m_Towers[itow].setRawAdc(value,idet);
m_Towers[itow].setEnergy( energy, idet );
//
// Accumulate sector-based statistics for this
// tower/detector ( 0=T, 1=P, 2=Q, 3=R ).
//
m_SectorStats[isec] -> addHit ( &m_Towers[itow], idet );
//
// Add to the list of hit towers and (if we exceed
// a user-specified threshold) the list of seed
// towers as well. Note that "HitTowers" means
// that the tower has responded with ADC > N sigma
// above threshold.
//
if ( idet == 0 && energy > 0. ) {
m_HitTowers.push_back( &m_Towers[itow] );
if ( energy >= m_SeedEnergy ) {
// Add to list of seed towers
m_SeedTowers.push_back( &m_Towers[itow] );
// Let this tower's neighbors know that it
// is a seed tower.
m_Towers[itow].setSeed();
}
}
} // !isSMD()
//--
//--
else if ( x -> isSMD() ) {
Int_t sector = x -> sec - 1;
Int_t plane = (Int_t)x -> plane - (Int_t)'U'; // 0=U, 1=V
Int_t strip = x -> strip - 1;
Float_t threshold = x -> thr;
//
// ADC must exceed the threshold specified w/in
// the database. If not, punt.
//
if ( value < threshold ) continue;
#ifndef ESMD_IDEAL_GAINS
if ( x -> gain <= 0. ) continue;
Float_t gain = x -> gain;
Float_t ped = x -> ped;
if ( gain <= 0 ) continue;
#endif
#ifdef ESMD_IDEAL_GAINS
Float_t gain = 23000.; // See StEEmcFastMaker
//Float_t gain = ideal_smd_gains[strip];
Float_t ped = x -> ped;
if ( gain <= 0 ) continue;
#endif
//
// If this strip's adc response is higher than the
// current high strip, set a new one...
//
if ( plane==0 ) {
if ( m_HighUStrip -> getRawAdc() > value )
m_HighUStrip = &m_Strips[ sector ] [ 0 ] [ strip ];
}
else {
if ( m_HighVStrip -> getRawAdc() > value )
m_HighVStrip = &m_Strips[ sector ] [ 1 ] [ strip ];
}
//
// 03/10/04 We're loading in gains by hand, which are expressed
// in terms of # MIPS / adc channel.
//
Float_t nmips = ( value - ped ) / gain; // number of mips
//$$$Float_t npe = ( value - ped ) / 8.0; // 1 p.e. destined for ch 8
// Number of photoelectrons / mip * number of mips
Float_t npe = ( 4.18 + 0.0107 * strip ) * nmips;
//
// For now, we do the following. Estimage dE/dx = 1.75 MeV/cm
// in the SMD strip. Height of the strip is 0.7 cm. However,
// an effective height of 0.55 cm in the "universal curve"
// should likely be assumed, since cuts were loose to get
// statistics generating that curve.
//
Float_t energy = nmips * 0.9625;
m_Strips[ sector ][ plane ][ strip ].setAdc( value - ped );
m_Strips[ sector ][ plane ][ strip ].setRawAdc( value );
m_Strips[ sector ][ plane ][ strip ].setEnergy( energy );
m_Strips[ sector ][ plane ][ strip ].setNumMips( nmips );
m_Strips[ sector ][ plane ][ strip ].setNumPhotoElectrons( npe );
//
// Accumulate sector statistics for this strip
//
m_SectorStats[sector] -> addHit ( &m_Strips[sector][plane][strip], plane );
if ( plane == 0 )
m_HitUStrips . push_back( &m_Strips[sector][plane][strip] );
else if ( plane == 1 )
m_HitVStrips . push_back( &m_Strips[sector][plane][strip] );
else
assert(0); // only 2 planes for now
}
} // loop over channels
} // loop over blocks
} // Fill()
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
// Fill using fzd files ( fzd --> EEeventDst )
#define SCALE_FACTOR 0.04
void EEezAnalysis::Fill(EEeventDst *event)
{
// As with the Fill(EEfeeRawEvent *event) method, except that it
// takes as input an EEeventDst, useful for direct processing
// of fzd files. Currently this method populates the tower
// channels only.
Int_t nSec = event->Sec->GetEntries();
for (Int_t isec=0; isec<nSec; isec++ ) {
EEsectorDst *secH =( EEsectorDst*)event->Sec->At(isec);
///////////////////////////////////////////////////////
//
// Fill towers
Int_t sec=secH->getID();
TClonesArray *hitA=secH->getTwHits();
assert(hitA);
// Loop over all hits
for ( Int_t ih=0; ih < hitA->GetEntries(); ih++ ) {
Float_t enerRaw;
Int_t eta;
Char_t sub;
EEtwHitDst *hit=(EEtwHitDst*)hitA->At(ih);
hit->get(sub,eta,enerRaw);
Float_t recoEnergy = enerRaw / SCALE_FACTOR;
Int_t isub = (Int_t)(sub-'A');
Int_t ieta = eta;
//<<<<<
Int_t itow = 60 * isec + 12 * isub + ieta;
assert( 0<= itow);
assert( itow<720);
m_Towers[itow].setAdc( (Float_t)(Int_t)( recoEnergy * ideal_gains[ieta] ), 0 );
m_Towers[itow].setEnergy( recoEnergy, 0 );
if ( recoEnergy > m_SeedEnergy ) {
m_SeedTowers.push_back( &m_Towers[itow] );
m_Towers[itow].setSeed();
}
}
//
///////////////////////////////////////////////////////
TClonesArray *hitPre1=secH->getPre1Hits();
assert(hitPre1);
TClonesArray *hitPre2=secH->getPre2Hits();
assert(hitPre2);
TClonesArray *hitPost=secH->getPostHits();
assert(hitPost);
///////////////////////////////////////////////////////
//
TClonesArray *hitSmdU=secH->getSmdUHits();
assert(hitSmdU);
for ( Int_t ih = 0; ih < hitSmdU -> GetEntries(); ih++ ) {
// Get the hit
EEsmdHitDst *hit = (EEsmdHitDst *)hitSmdU -> At(ih);
Int_t strip = hit -> strip() - 1;
Float_t eraw = hit -> energy();
m_Strips[sec-1][0][strip].setEnergy(eraw/1.975);
m_Strips[sec-1][0][strip].setAdc(eraw);
}
///////////////////////////////////////////////////////
//
TClonesArray *hitSmdV=secH->getSmdVHits();
assert(hitSmdV);
for ( Int_t ih = 0; ih < hitSmdV -> GetEntries(); ih++ ) {
// Get the hit
EEsmdHitDst *hit = (EEsmdHitDst *)hitSmdV -> At(ih);
Int_t strip = hit -> strip() - 1;
Float_t eraw = hit -> energy();
m_Strips[sec-1][1][strip].setEnergy(eraw/1.975);
m_Strips[sec-1][1][strip].setAdc(eraw);
}
}
}
#endif
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
TH2F EEezAnalysis::getHistogram()
{
// Creates a 2D eta vs phi plot from the tower readouts.
TH2F hist = TH2F("eeEzAnaly","Energy vs eta vs phi",60,0.,60.,12,0.,12.);
for ( Int_t i = 0; i < 720; i++ ) {
Int_t eta = i % 12;
Int_t phi = i / 12;
hist.SetBinContent(phi,eta, m_Towers[i].getEnergy() );
}
return hist;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::Init()
{
// Main initialization routine.
Info("Init","Initializing EEezAnalysis");
//$$$ InitStrips();
//$$$ InitTowers();
//$$$ InitTowerStrips();
InitGeometry();
// Start with high-tower and high-strip pointers
// pointing someplace, doesn't matter where.
m_HighTower = &m_Towers[0];
m_HighUStrip = &m_Strips[0][0][0];
m_HighVStrip = &m_Strips[0][1][0];
#ifndef __ROOT__
for ( Int_t i = 0; i < 12; i++ ) {
m_SectorStats[i] -> setSeedEnergy( m_SeedEnergy );
}
#endif
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::InitGeometry()
{
//
// Initialize the geometrical relations between towers and strips
//
InitStrips();
InitTowers();
InitTowerStrips();
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::InitTowerStrips()
{
//
// For each tower, access the EEmcSmdMap to determine the range
// of U and V strips which are nominally associated with each
// tower. It then calls addUStrip() and addVStrip() on the tower
// to add pointers to these strips.
//
// Additionally, access the EEmcSmdMap to determine the list of
// towers which share the U and V strips. These towers are
// added to the SMD-shared-tower lists kept by each tower.
//
for ( Int_t i = 0; i < 720; i ++ ) {
// Get the sector, subsector and etabin for this tower
Int_t sector = m_Towers[i] . getSector();
Int_t subsector = m_Towers[i] . getSubSector();
Int_t etabin = m_Towers[i] . getEtabin();
// Get the min and max U- and V-strip indices for this tower
Int_t umin=-1,umax=-1;
Int_t vmin=-1,vmax=-1;
m_EEmcSmdMap -> getRangeU( sector, subsector, etabin, umin, umax );
m_EEmcSmdMap -> getRangeV( sector, subsector, etabin, vmin, vmax );
// Now, loop over all strips within the specified ranges,
// and add them to the list of strips within the fiducial
// volume of this tower.
for ( Int_t u = umin; u <= umax; u++ ) {
m_Towers[i] . addUStrip ( &m_Strips[sector][0][u] );
}
for ( Int_t v = vmin; v <= vmax; v++ ) {
m_Towers[i] . addVStrip ( &m_Strips[sector][1][v] );
}
}
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::InitStrips()
{
// Initializes the m_Strips array
for ( Int_t i = 0; i < 12; i++ ) {
for ( Int_t j = 0; j < 288; j++ ) {
for ( Int_t k = 0; k < 2; k++ ) {
// Note the odd ordering of ikj... for sector plane strip
m_Strips[i][k][j].clear();
m_Strips[i][k][j].setIndex(j);
m_Strips[i][k][j].setSector(i);
m_Strips[i][k][j].setPlane(k);
}
}
}
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::RestoreIndices()
{
// Reinitializes the tower indices if they get lost. Shouldn't be needed.
// Initialize tower index
for ( Int_t itow = 0; itow < 720; itow++ )
m_Towers[itow].setIndex(itow);
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::InitTowers()
{
// Initializes the geometric description of the towers in the m_Towers
// array. Loops over all towers, and sets their index, and sets pointers
// to all neighboring towers.
// Loop over all towers and initialize pointers
// to neighboring towers
// Clear all towers and initialize tower index
for ( Int_t itow = 0; itow < 720; itow++ ) {
m_Towers[itow].clear();
m_Towers[itow].setIndex(itow);
}
// loop over all bins in phi
// (code below was originally written to deal with m_Towers[12][5][12]...
// cut and paste here...)
for ( Int_t iphi = 0; iphi < 60; iphi++ ) {
Int_t isec = iphi / 5;
Int_t isub = iphi % 5;
// loop over all eta bins
for ( Int_t ieta = 0; ieta < 12; ieta++ ) {
// Loop over adjacent phi bins
for ( Int_t jphi = iphi - 1; jphi <= iphi + 1; jphi++ ) {
Int_t jsec = (jphi / 5);
Int_t jsub = (jphi % 5);
// Remap sector and subsector at boundaries
if ( jsub == -1 ) {
jsub = 4;
jsec = 11;
}
if ( jsec == 12 ) {
jsec = 0;
jsub = 0;
}
// Loop over adjacent eta bins
for ( Int_t jeta = ieta - 1; jeta <= ieta + 1; jeta++ ) {
// Avoid self
if ( isec == jsec && isub == jsub && ieta == jeta ) continue;
// Add pointer to neighboring tower, or set NULL if
// out of the eta boundary
if ( jeta >= 0 && jeta < 12 ) {
m_Towers[60*isec + 12*isub + ieta].addNeighbor( &m_Towers[60*jsec + 12*jsub +jeta] );
}
} // adjacent eta
} // adjacent phi
//$$$ m_Towers[60*isec + 12*isub + ieta].printNeighbors();
} // eta
} // phi
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::Save ( EEezAnalysis *analysis )
{
// Save() provides a method for copying the current state of the
// EEezAnalysis object. It may or may not be useful at some point,
// and is kept for now.
// Given a pointer to another instance of the analysis class,
// save all of the struck towers within this class. This
// would be useful for storing an event for background mixing.
// In particular, we need to be certain that the data in the
// tower classes gets copied over without getting the pointers
// to adjacent towers wrong.
// Clear this instance from any previous analysis
Clear();
// Get the seed energy from the other analysis
m_SeedEnergy = analysis -> getSeedEnergy();
//
// Copy the hit towers to local storarge.
//
EEezTowerPtrVecIter_t iter = analysis -> getHitTowers() -> begin();
while ( iter != analysis -> getHitTowers() -> end() ) {
Int_t index = (*iter) -> getIndex();
// Copy over energy for preshower, postshower and tower
for ( Int_t i = 0; i < 3; i++ ) {
m_Towers[ index ].setEnergy( (*iter) -> getEnergy(i), i );
m_HitTowers.push_back ( &m_Towers[ index ] );
}
// If this tower satisfies the seed energy, add it to the
// list of seed towers
if ( (*iter) -> getEnergy(0) > m_SeedEnergy ) {
m_SeedTowers.push_back( &m_Towers[ index ] );
}
iter++;
}
//
// NOTE: Copy over strips, when they become available...
//
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::print()
{
// Will do something eventually.
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::addTowerHit ( Int_t sec, Int_t sub, Int_t eta, Float_t adc, Int_t idet )
{
// Public method used to add a tower hit. Handles ADC --> energy
// conversion, storing seed towers, etc... It's a public method so
// that external wrapper functions may be used to fill the m_Towers
// array for each event.
//
// Will eventually be moved into the loop in "Fill".
assert(sec<12); //-- Indices start from 0 in c++
// Index of the tower in the array
Int_t itow = 60 * sec + 12 * sub + eta;
assert(itow<720); //-- We're just out of bounds today.
// Get the pedestal and gain for this tower (pre/post) from
// the database
static Char_t subsectors[] = { 'A','B','C','D','E' };
static Char_t detectors[] = { 'T', 'P', 'Q', 'R' };
#if 0
std::cout << ">>>>>>>>>>>";
std::cout << "Getting tail (tail?) -- " << sec+1
<< " " << subsectors[sub] << " "
<< eta + 1
<< std::endl << std::flush;
#endif
// Get the database entry for this detector (note, DB expects
// sectors, eta counted from 1 not 0 for conformity with Star's ...
// conventions.)
const EEmcDbItem_t *dbitem = m_EEmcDb -> getTile( sec+1,subsectors[sub],eta+1, detectors[idet] );
assert(dbitem);
// Check if null
if ( dbitem == 0 ) return;
// Check that the tower (pre/post) is valid
if ( dbitem -> fail ) return;
// Mask off user-requested status bits
if ( m_StatMask & dbitem -> stat ) return;
// Correct ADC for pedestal and gain
Float_t value = adc;
Float_t ped = dbitem -> ped;
Float_t gain = dbitem -> gain;
Float_t threshold = dbitem -> thr;
// Verify that the ADC response is above the threshold
// specified w/in the database
if ( value < threshold ) return;
//--
//-- NOTE: Make no provision for "ideal" gains for the
//-- time being. Just trap gain == 0.
//--
if ( gain <= 0. ) return;
Float_t energy = ( value - ped ) / gain;
m_Towers[itow].setAdc( value-ped, idet );
m_Towers[itow].setRawAdc(value,idet);
m_Towers[itow].setEnergy( energy, idet );
#ifndef __ROOT__
// Accumulate sector-based statistics for this
// tower/detector ( 0=T, 1=P, 2=Q, 3=R ).
m_SectorStats[sec] -> addHit ( &m_Towers[itow], idet );
#endif
// Add to the list of hit towers and (if we exceed
// a user-specified threshold) the list of seed
// towers as well. Note that "HitTowers" means
// that the tower has responded with ADC > N sigma
// above threshold.
if ( idet == 0 && energy > 0. ) {
m_HitTowers.push_back( &m_Towers[itow] );
if ( energy >= m_SeedEnergy ) {
// Add to list of seed towers
m_SeedTowers.push_back( &m_Towers[itow] );
// Let this tower's neighbors know that it
// is a seed tower.
m_Towers[itow].setSeed();
}
}
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezAnalysis::addSmdHit ( Int_t sec, Int_t plane, Int_t strip, Float_t adc )
{
// Adds a hit to the specified SMD strip. Handles ADC --> energy
// conversion, storing seed towers, etc... It's a public method so
// that external wrapper functions may be used to fill the m_Towers
// array for each event.
//
// Will eventually be moved into the loop in "Fill".
assert(sec<12); // c++ indices are from 0!
assert(plane<2); // U = 0, V = 1
assert(strip<288); // strips numbered from 0-287
static Char_t planes[] = { 'U', 'V' };
// Get the pedestal and gain for this SMD strip from
// the database
// Access the database. Note that the DB counts from 1, not zero.
const EEmcDbItem_t *dbitem =
m_EEmcDb -> getByStrip ( sec+1, planes[plane], strip+1 );
// Check if null
if ( dbitem == 0 ) return;
// Check that the strip is valid
if ( dbitem -> fail ) return;
// Mask off user-requested status bits
if ( m_StatMask & dbitem -> stat ) return;
// Correct ADC for pedestal and gain
Float_t value = adc;
Float_t ped = dbitem -> ped;
Float_t gain = dbitem -> gain;
Float_t threshold = dbitem -> thr;
#if 0
gain = 23000 * ( 0.9625 / 1000.0 );
threshold = 0.000;
ped = 0.000;
#endif
// Verify that the ADC response is above the threshold
// specified w/in the database
if ( value < threshold ) return;
//--
//-- NOTE: Make no provision for "ideal" gains for the
//-- time being. Just trap gain == 0.
//--
if ( gain <= 0. ) return;
Float_t nmips = ( value - ped ) / gain; // number of mips
Float_t npe = ( 4.18 + 0.0107 * strip ) * nmips;
#if 0
if ( nmips < 1.0 ) return;
#endif
// For now, we do the following. Estimage dE/dx = 1.75 MeV/cm
// in the SMD strip. Height of the strip is 0.7 cm. However,
// an effective height of 0.55 cm in the "universal curve"
// should likely be assumed, since cuts were loose to get
// statistics generating that curve.
Float_t energy = nmips * 0.9625;
m_Strips[ sec ][ plane ][ strip ].setAdc( value - ped );
m_Strips[ sec ][ plane ][ strip ].setRawAdc( value );
m_Strips[ sec ][ plane ][ strip ].setEnergy( energy );
m_Strips[ sec ][ plane ][ strip ].setNumMips( nmips );
m_Strips[ sec ][ plane ][ strip ].setNumPhotoElectrons( npe );
#ifndef __ROOT__
// Accumulate sec statistics for this strip
m_SectorStats[ sec ] -> addHit ( &m_Strips[ sec ][plane][strip], plane );
#endif
if ( plane == 0 )
m_HitUStrips . push_back( &m_Strips[ sec ][plane][strip] );
else if ( plane == 1 )
m_HitVStrips . push_back( &m_Strips[ sec ][plane][strip] );
else
assert(0); // only 2 planes for now
}
void EEezAnalysis::setDb( EEmcDb_t *db )
{
// Initialize the pointer to the eemc database maker
assert(db);
m_EEmcDb = db;
}
/////////////////////////////////////////////////////////////////////////////
// COMPILE ONLY IN STROOT ENVIRONMENT ///////////////////////////////////////
#if __ROOT__
Int_t EEezAnalysis::Make( StMuEmcCollection *emc )
{
// Loop over all hits in the EMC collection and add them to our
// internal array of EEezTowers.
assert(emc); //-- Need the data after all
//-- Loop over all 720 towers
for ( Int_t ihit = 0; ihit < emc -> getNEndcapTowerADC(); ihit++ ) {
Int_t adc, sec, sub, eta;
emc -> getEndcapTowerADC ( ihit, adc, sec, sub, eta );
sec--; // Counting from 1 is insane in c++! It ends here.
sub--; // Eveerything in my classes assume indexes from 0.
eta--; //
assert( sec >= 0 && sec < 12 ); // Indexing errors detected
assert( sub >= 0 && sub < 5 ); // Indexing errors detected
assert( eta >= 0 && eta < 12 ); // Indexing errors detected
addTowerHit(sec,sub,eta,adc,0);
}
//-- Loop over all preshower/postshower hits
for ( Int_t ihit = 0; ihit < emc -> getNEndcapPrsHits(); ihit++ ) {
Int_t adc, sec, sub, eta, det;
StMuEmcHit *hit = emc -> getEndcapPrsHit(ihit, sec, sub, eta, det);
adc = hit -> getAdc();
sec--; // Counting from 1 is insane in c++! It ends here.
sub--; // Eveerything in my classes assume indexes from 0.
eta--; //
assert( sec >= 0 && sec < 12 ); // Indexing errors detected
assert( sub >= 0 && sub < 5 ); // Indexing errors detected
assert( eta >= 0 && eta < 12 ); // Indexing errors detected
assert( det >= 1 && det < 4 ); // Indexing errors detected
addTowerHit(sec,sub,eta,det);
}
//-- Loop over all SMD hits
Char_t cpl[] = { 'U','V' };
for ( Int_t plane = 0; plane < 2; plane++ )
for ( Int_t ihit = 0; ihit < emc -> getNEndcapSmdHits(cpl[plane]); ihit++ ) {
Int_t sec, strip, adc;
StMuEmcHit *hit = emc -> getEndcapSmdHit( cpl[plane], ihit, sec, strip );
adc = hit -> getAdc();
sec--; // Counting from 1 is insane in c++! It ends here.
strip--; // Eveerything in my classes assume indexes from 0.
assert( sec >= 0 && sec < 12 ); // Indexing errors detected
assert( strip >= 0 && strip < 288 ); // Indexing errors detected
addSmdHit(sec,plane,strip,adc);
}
return kStOK;
}
#endif
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.