/////////////////////////////////////////////////////////////////////////////
//
// EEsmdCalibration
//
// Author: Jason C. Webb <jwebb@iucf.indiana.edu>
//
// This class is designed to iteratively measure the ESMD strip-to-strip
// gains....
//
// Ungh... It's looking more and more like a fortran program.
//
/////////////////////////////////////////////////////////////////////////////
#include "EEsmdCalibration.h"
#include "EEezAnalysis.h"
#include "EEezClAnalysis.h"
#include "EEmcDb.h"
#include "EEsmdProfile.h"
#include "EEezCluster.h"
#include "StEEmcUtil/EEfeeRaw/EEname2Index.h"
#include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
#include "StEEmcDbMaker/EEmcDbItem.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TMath.h"
#include "TVector3.h"
#include <iostream>
ClassImp(EEsmdCalibration);
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
EEsmdCalibration::EEsmdCalibration( const Char_t *name, const Char_t *title )
: TDirectory(name,title)
{
// Class constructor
m_NSigma = 3.0;
m_YieldMin = 15.0;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
Int_t EEsmdCalibration::Init()
{
// Initialization of histograms and internal data structures.
assert(m_EEezAnaly); // need pointer to the analysis class
assert(m_EEmcDbase); // need pointer to database
assert(m_EEezClust); // need pointer to cluster finder
assert(m_EEsmdProf); // need pointer to smd profiler
// Book histograms to hold the mean gain shifts required to "smooth"
// out the SMD U and V distributions
cd();
m_NumClusters = new TH1F("m_NumClusters","Total number of tower-clusters accumulated in each sector",12,0.,12.);
// Loop over all sectors/both planes
for ( Int_t i = 0; i < 12; i++ ) {
TString sec = "";
sec += ( i<9 ) ? "0" : "";
sec += (i+1);
for ( Int_t iuv = 0; iuv < 2; iuv++ ) {
TString plane = sec;
plane += (iuv) ? "V" : "U";
TString name = "gshift"; name += plane;
TString title = "Mean, fractional residuals for "; title += plane;
m_GainShift[i][iuv] = new TProfile(name,title,288,0.,288.,-2.,2.);
name = "sumMip"; name += plane;
title = "Total number of mips (after cuts) vs strip ID";
m_SumMips[i][iuv] = new TH1F(name,title,288,0.,288.);
name = "sumHit"; name += plane;
title = "Total number of hits (after cuts) vs strip ID";
m_SumHits[i][iuv] = new TH1F(name,title,288,0.,288.);
name = "crossCalib"; name += plane;
title = "Ratio of N_{mip} for max strip in plane to N_{mip}^{max} in other plane";
m_CrossCalib[i][iuv] = new TProfile(name,title,288,0.,288.,0.,100.);
name = "meanMip"; name += plane;
title = "Mean number of mips vs strip index";
m_MeanMips[i][iuv] = new TProfile(name,title,288,0.,288.,0.,1000.);
}
}
m_Chi2VsStrip = new TH2F("m_Chi2VsStrip","#chi^2 / ndf vs strip index of strips w/in 3#sigma of mean",288,0.,288.,50,0.,10.);
m_Chi2VsEclust = new TH2F("m_Chi2VsEclust","#chi^2 / ndf vs E_{c}",30,1.0,4.0,50,0.,10.);
m_YieldVsEclust = new TProfile("m_YieldVsEclust","N_{mip} vs E_{c}",30,1.0,4.0,0,1000.);
m_SigmaVsEclust = new TProfile("m_SigmaVsEclust","Width of fit vs E_{c}",30,1.,4.,0.,10.);
m_FracVsEclust = new TProfile("m_FracVsEclust","Fraction of 2nd gaussian vs E_{c}",30,1.,4.,0.,1.);
m_Chi2VsShape = new TH2F("m_Chi2VsShape","#chi^2 / ndf vs eSeed/eCluster",30,0.7,1.0,50,0.,10.);
m_CutEvents = new TH1F("m_CutEvents","Cut event distribution",1,0,1);
m_CutEvents -> SetBit( TH1::kCanRebin );
//
// Residual spectra for specific strips
//
Int_t
strips[] = { 24, 25, 26,
49, 50, 51,
74, 75, 76,
149, 150, 151 };
for ( Int_t i = 0; i < 288; i++ ) m_SelectStripMap[i] = -1;
for ( Int_t i = 0; i < 12; i++ ) {
m_SelectStrips[i] = strips[i];
TString name = "strip"; name += strips[i];name += "resid";
//$$$ --> Moverd = new TH1F(name,"Selected strip residual spectrum",200,-4.,4.);
name += "VsEclust";
m_ResidualsVsEclust[i] = new TProfile(name,"Selected strip residuals vs E_{cluster}",30,1.,4.,-10.,10.);
name.ReplaceAll("Eclust","Chi2");
m_ResidualsVsChi2[i] = new TProfile(name,"Selected strip residuals vs #chi^2 of fit",30,0.,10.,-10.,10.);
m_SelectStripMap[ strips[i]-1 ] = i;
name.ReplaceAll("Chi2","Dist");
m_ResidualsVsDist[i] = new TProfile(name,"Selected strip residuals vs distance from mean (num strips)",40,-5.,5.,-10.,10.);
}
// Store residuals for every strip in sub dir
TDirectory *sub1 = new TDirectory("residuals","Residuals of individual strips");
sub1 -> cd();
for ( Int_t sec = 0; sec < 12; sec++ ) {
for ( Int_t str = 0; str < 288; str++ ) {
TString name = "m_Residuals";
name += "["; name += sec; name += "]";
name += "["; name += str; name += "]";
m_Residuals[sec][str] = new TH1F(name,"Fractional strip residuals ( U -0.5, V ++0.5 )",50,-1.,1.);
}
}
sub1 -> cd("..");
// Store yield to fit ratios in subdirectory
TDirectory *sub = new TDirectory("yield2fit","Yield-to-fit ratios");
sub -> cd();
// Loop over all sectors/both planes
for ( Int_t i = 0; i < 12; i++ ) {
for ( Int_t iuv = 0; iuv < 2; iuv++ ) {
for ( Int_t strip = 0; strip < 288; strip++ ) {
TString sec;
sec += "["; sec += i; sec += "]";
sec += "["; sec += iuv; sec += "]";
sec += "["; sec += strip; sec += "]";
TString name = "m_YieldToFit"; name += sec;
TString title = "Yield-to-fit ratio";
m_YieldToFit[i][iuv][strip] = new TH1F(name,title,100,0.,2.);
m_YieldToFit[i][iuv][strip] -> Sumw2();
}
}
}
sub -> cd("..");
// Return to previous directory
cd("..");
return 1;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
Int_t EEsmdCalibration::Make()
{
// Process one event. NOTE... given root's terse documentation
// of TH1::FindBin(), there is a nonzero chance of a binning error
// here... if rumors of underflow and overflow bins turn out to
// be true.
cd();
// XY distribution of clusters with no U or no V strips under seed
static TH2F *h_make_xy_nou;
static TH2F *h_make_xy_nov;
static TH1F *h_make_pre1_adc;
static TH1F *h_make_pre2_adc;
static TH1F *h_make_post_adc;
static TH1F *h_make_post_nosmd_adc;
static Int_t first = 1;
if ( first ) {
first = 0;
TDirectory *sub = new TDirectory("makeCuts","XY distributions of events eliminated in Make()");
sub -> cd();
h_make_xy_nou = new TH2F("h_make_xy_nou","XY distribution of clusters with no seed-tower SMD-U response",100,-175.,175.,100,-175.,175.);
h_make_xy_nov = new TH2F("h_make_xy_nov","XY distribution of clusters with no seed-tower SMD-U response",100,-175.,175.,100,-175.,175.);
h_make_pre1_adc = new TH1F("h_make_pre1_adc","Preshower 1--ADC spectrum for isolated seed towers",1024,0.,1024.);
h_make_pre2_adc = new TH1F("h_make_pre2_adc","Preshower 2--ADC spectrum for isolated seed towers",1024,0.,1024.);
h_make_post_adc = new TH1F("h_make_post_adc","Postshower--ADC spectrum for isolated seed towers",1024,0.,1024.);
h_make_post_nosmd_adc = new TH1F("h_make_post_nosmd_adc","Postshower--ADC spectrum for isolated seed towers w/ no SMD response",1024,0.,1024.);
sub -> cd("..");
}
//-- Loop over all clusters
for ( Int_t ic = 0; ic < m_EEezClust -> getNClusters(); ic++ ) {
// Get a pointer to the cluster
EEezCluster *cluster = m_EEezClust -> getClusterPtr(ic);
m_NumClusters -> Fill ( (Float_t)cluster -> getSeedTower() -> getSector() + 0.5 );
// Get a pointer to the seed tower as well
EEezTower *seed = cluster -> getSeedTower();
h_make_pre1_adc -> Fill ( seed -> getAdc(1) );
h_make_pre2_adc -> Fill ( seed -> getAdc(2) );
h_make_post_adc -> Fill ( seed -> getAdc(3) );
// Keep track of the total number of identified clusters
m_CutEvents -> Fill ( "N_CLUSTERS", 1.0 );
// Check if it's isolated, and continue if not
// if ( isolationCut( cluster ) ) {
if ( isolationCut( ic ) ) {
m_CutEvents -> Fill ( "FAIL_ISOLATION_CUT", 1.0 );
continue;
}
// Cumulate total number of isolated clusters
m_CutEvents -> Fill( "N_ISO_CLUSTERS", 1.0 );
// Get the energy of the cluster for future use
Float_t eCluster = cluster -> getEnergy();
Float_t eSeed = cluster -> getSeedTower() -> getEnergy();
// Get the sector we are dealing with
Int_t sector = cluster -> getSeedTower() -> getSector();
// Get the SMD U and V distributions
TH1F *udist = m_EEsmdProf -> getProfileU(ic);
TH1F *vdist = m_EEsmdProf -> getProfileV(ic);
//
// Punt if the histogram has too few entries (code crashes
// w/out this cut!).
//
TVector3 momentum = cluster -> getMomentum();
TVector3 direction = momentum.Unit();
Int_t fail = 0;
if ( udist -> GetEntries() < 1 ) {
m_CutEvents -> Fill ( "FAIL_U_DIST_EMPTY", 1.0 );
h_make_xy_nou ->Fill( (kEEmcZSMD*direction).X(), (kEEmcZSMD*direction).Y() );
fail++;
}
if ( vdist -> GetEntries() < 1 ) {
m_CutEvents -> Fill ( "FAIL_V_DIST_EMPTY", 1.0 );
h_make_xy_nov ->Fill( (kEEmcZSMD*direction).X(), (kEEmcZSMD*direction).Y() );
fail++;
}
if ( fail ) {
h_make_post_nosmd_adc -> Fill ( seed -> getAdc(3) );
continue;
}
fail = 0;
// Cumulate number of isolated clusters w/ U and V profiles
m_CutEvents -> Fill( "N_ISO_CLUSTERS_WITH_SMD", 1.0 );
// Retrieve the fits to the SMD U and V distributions
TF1 *ufit = (TF1 *)udist -> GetListOfFunctions() -> At(0);
TF1 *vfit = (TF1 *)vdist -> GetListOfFunctions() -> At(0);
// Get the fit-yields, means and widths of the two distributions
Float_t uyield = ufit -> GetParameter(0);
Float_t vyield = vfit -> GetParameter(0);
Float_t umean = ufit -> GetParameter(1);
Float_t vmean = vfit -> GetParameter(1);
Float_t usigma = ufit -> GetParameter(2);
Float_t vsigma = vfit -> GetParameter(2);
Float_t uchi2 = ufit -> GetChisquare();
Float_t vchi2 = vfit -> GetChisquare();
Float_t undf = ufit -> GetNDF();
Float_t vndf = vfit -> GetNDF();
Float_t ufrac = ufit -> GetParameter(3);
Float_t vfrac = vfit -> GetParameter(3);
// Fill chi2 vs ecluster and other parameters
if ( undf ) {
m_Chi2VsEclust -> Fill ( eCluster, uchi2 / undf );
m_Chi2VsShape -> Fill ( eSeed/eCluster, uchi2 / undf );
m_YieldVsEclust -> Fill ( eCluster, uyield );
m_SigmaVsEclust -> Fill ( eCluster, usigma );
m_FracVsEclust -> Fill ( eCluster, ufrac );
}
if ( vndf ) {
m_Chi2VsEclust -> Fill ( eCluster, vchi2 / vndf );
m_Chi2VsShape -> Fill ( eSeed/eCluster, vchi2 / vndf );
m_YieldVsEclust -> Fill ( eCluster, vyield );
m_SigmaVsEclust -> Fill ( eCluster, vsigma );
m_FracVsEclust -> Fill ( eCluster, vfrac );
}
// When we loop over the bins, we will not keep residuals
// for the bin containing the mean or the two strips
// adjacent to it. These strips provide the most constraints
// in the fit.
Int_t umeanbin = udist -> FindBin(umean);
Int_t vmeanbin = udist -> FindBin(vmean);
// We will loop over some number of bins, defined by how many
// sigma away from the mean we have requested...
Int_t umin = umeanbin - 4; // udist -> FindBin( umean - 4.5 );
Int_t umax = umeanbin + 4; // udist -> FindBin( umean + 4.5 );
Int_t vmin = vmeanbin - 4; // vdist -> FindBin( vmean - 4.5 );
Int_t vmax = vmeanbin + 4; // vdist -> FindBin( vmean + 4.5 );
// Now, calculate the fractional residual for each
// strip, and fill the profile histogram... NOTE:
// bin number here should equal strip number, which
// is indexed from 0...
// NOTE: Forget normalization for now..., renorm at the end
// Don't forget... umin and umax are bin numbers! Likely
// just 0 to nmax.... so we need to add the lower-edge
// of the first bin when we fill the histograms (and another
// 0.5 so we're in the center of the bin).
Float_t xmin = udist -> GetBinLowEdge(0);
Int_t umaxbin = -1;
Float_t nmip_max = 0.;
//--
//-- Before binning, verify that the event satisfies
//-- both simple cuts on the SMD profile histograms,
//-- and that it satisfies some QA cuts on the
//-- fits. The logic gets a little convoluted here.
//--
Int_t goProfileU = ( profileQAcut(ic,0) == 0 );
Int_t goProfileV = ( profileQAcut(ic,1) == 0 );
if ( !goProfileU ) {
m_CutEvents -> Fill ( "FAIL_PROFILE_QA_U", 1.0 );
}
else if ( !goProfileV ) {
m_CutEvents -> Fill ( "FAIL_PROFILE_QA_V", 1.0 );
}
if ( goProfileU && goProfileV ) {
m_CutEvents -> Fill ( "N_CLUSTERS_WITH_PROFILE_QA", 1.0 );
}
else
continue;
Int_t goFitU = ( fitQAcut(ic,0) == 0 );
Int_t goFitV = ( fitQAcut(ic,1) == 0 );
if ( !goFitU ) {
m_CutEvents -> Fill ( "FAIL_FIT_QA_U", 1.0 );
}
else if ( !goFitV ) {
m_CutEvents -> Fill ( "FAIL_FIT_QA_V", 1.0 );
}
if ( goFitU && goFitV ) {
m_CutEvents -> Fill ( "N_CLUSTERS_WITH_FIT_QA", 1.0 );
}
else
continue;
//
// Loop over the ~9 strips near the mean strip
//
for ( Int_t u = umin; u <= umax; u++ ) {
// Calculate x at the center of the strip
Float_t x = udist -> GetBinLowEdge(u) + 0.5 * udist -> GetBinWidth(u);
// Get the mip yield for this bin
Float_t d = udist -> GetBinContent(u);
// Calculate the prediction from the fit
Float_t f = ((TF1*)udist -> GetListOfFunctions() -> At(0)) -> Eval( x );
// Get the chi2 from the fit
// Determine the residual
Float_t r = d - f;
// And estimate the error on the residual
Float_t e = udist -> GetBinError(u);
if ( d > nmip_max ) { nmip_max = d; umaxbin = u; }
if ( !f ) continue;
Float_t frac = d / f;
// Fill the gain-shift histograms (TProfile of the mean-
// fractional residual versus strip index)
if ( d )
m_GainShift[sector][0] -> Fill(x,frac);
if ( undf * d )
m_Chi2VsStrip -> Fill ( x, uchi2 / undf );
if ( f*d>0. ) {
// Fill with the ratio, weighted by 1/sigma**2
// m_YieldToFit[sector][0][u-1] -> Fill( d/f, f*f/d );
m_YieldToFit[sector][0][u-1] -> Fill( d/f );
}
m_SumMips[sector][0] -> Fill(x,d);
m_SumHits[sector][0] -> Fill(x);
m_MeanMips[sector][0] -> Fill(x,d);
m_Residuals[sector][u] -> Fill(x,frac-0.5);
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
xmin = vdist -> GetBinLowEdge(0);
Int_t vmaxbin = -1;
nmip_max = 0.;
for ( Int_t v = vmin; v <= vmax; v++ ) {
//$$$Float_t x = (Float_t)v+xmin+0.5;
Float_t x = vdist -> GetBinLowEdge(v) + 0.5 * vdist -> GetBinWidth(v);
//Float_t r = vresid -> GetBinContent(v);
Float_t d = vdist -> GetBinContent(v);
Float_t f = ((TF1*)vdist -> GetListOfFunctions() -> At(0)) -> Eval( x );
Float_t r = d - f;
//Float_t e = vresid -> GetBinError(v);
Float_t e = vdist -> GetBinError(v);
if ( d > nmip_max ) {
nmip_max = d;
vmaxbin = v;
}
if ( !f ) continue;
Float_t frac = r/f;
if ( d ) {
m_GainShift[sector][1] -> Fill( x, frac );
}
if ( f*d>0. ) {
// Fill with the ratio, weighted by 1/sigma**2
//$$$ m_YieldToFit[sector][1][v-1] -> Fill( d/f, f*f/d );
m_YieldToFit[sector][1][v-1] -> Fill ( d/f );
}
if ( vndf * d ) m_Chi2VsStrip -> Fill ( x, vchi2 / vndf );
m_SumMips[sector][1] -> Fill(x,d);
m_SumHits[sector][1] -> Fill(x);
m_MeanMips[sector][1] -> Fill(x,d);
m_Residuals[sector][v] -> Fill(x,frac+0.5);
}
if ( umin <= umaxbin && vmin <= vmaxbin &&
umax >= umaxbin && vmax >= vmaxbin ) {
Float_t umip = udist -> GetBinContent(umaxbin);
Float_t vmip = udist -> GetBinContent(vmaxbin);
Float_t ux = udist -> GetBinLowEdge(0) + (Float_t)umaxbin + 0.5;
Float_t vx = vdist -> GetBinLowEdge(0) + (Float_t)vmaxbin + 0.5;
m_CrossCalib[sector][0] -> Fill(ux,umip/vmip);
m_CrossCalib[sector][1] -> Fill(vx,vmip/umip);
}
}
cd("..");
return 1;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEsmdCalibration::Clear()
{
// Clear internal data structures
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEsmdCalibration::newDb( Int_t ifirst, Int_t ilast )
{
// Calculates and loads into the database structures a new set of
// gains, based on the current state of the profile histograms.
assert(m_EEmcDbase);
Info("newDb","Setting database values from profile histograms");
for ( Int_t i = ifirst; i <= ilast; i++ ) {
TString sec = "";
if ( (i+1)<10 ) sec+= "0";
sec += (i+1);
for ( Int_t j = 0; j < 2; j++ ) {
TString plane = sec;
plane += ( j ) ? "V" : "U";
TString histo = "gshift"; histo += plane;
cd();
TH1F *h = (TH1F *) Get(histo);
cd("..");
std::cout << "Processing plane = " << plane << std::endl << std::flush;
Float_t oldSum = 0.;
Float_t newSum = 0.;
for ( Int_t k = 0; k < 288; k++ ) {
TString name = plane;
if ( (k+1)<100 ) name += "0";
if ( (k+1)<10 ) name += "0";
name += (k+1);
std::cout << name << " " << std::flush;
Int_t index = EEname2Index( name );
std::cout << index << " " << std::flush;
const EEmcDbItem *x = m_EEmcDbase -> getByIndex( index );
if ( !x -> fail )
oldSum += x -> gain; // initial gain sum for this plane
std::cout << " gain = " << x -> gain << " --> " << std::flush;
// ---<<< buggy software >>>---
//$$$ Float_t shift = h -> GetBinContent(k);
//$$$ Float_t gain = x->gain / ( 1 + shift );
Float_t shift = h -> GetBinContent( k+1 ); // bins numbered from 1
Float_t gain = x->gain * ( 1.0 + shift ); // pos. residual = low gain
std::cout << gain << std::endl << std::flush;
// Direct, naughty access to the database gain
if ( TMath::Abs(shift) < m_MaxShift )
(m_EEmcDbase -> byIndex + index) -> gain = gain;
if ( !x -> fail )
newSum += x -> gain; // new gain sum for this plane
}
// Preserve normalization of the planes
#define PRESERVE_NORMALIZATION
#ifdef PRESERVE_NORMALIZATION
for ( Int_t k = 0; k < 288; k++ ) {
TString name = plane;
if ( (k+1)<100 ) name += "0";
if ( (k+1)<10 ) name += "0";
name += (k+1);
Int_t index = EEname2Index( name );
const EEmcDbItem *x = m_EEmcDbase -> getByIndex( index );
if ( !x -> fail ) {
Float_t gain = x->gain * oldSum / newSum;
(m_EEmcDbase -> byIndex + index) -> gain = gain;
}
}
#endif
}
}
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEsmdCalibration::kludgeGains()
{
// Sets initial SMD gains IN THE DATABASE to a parameterization of the
// mean #pe per MIP, setting the single photoelectron peak at ADC
// channel 8.
for ( Int_t i = 1; i <= 12; i++ ) {
TString sec = "";
if ( i < 10 ) sec += "0";
sec += i;
for ( Int_t j = 0; j < 2; j++ ) {
TString plane = sec;
plane += (j) ? "V" : "U";
for ( Int_t k = 1; k <= 288; k++ ) {
TString name=plane;
if ( k < 100 ) name += "0";
if ( k < 10 ) name += "0";
name += k;
Int_t index = EEname2Index(name);
//
// ----<<<<< temporary kludge for gains >>>>>----
Int_t stripID = k;
Float_t y=0;
Float_t ya=0,yb=0,xa=1,xb=2;
if ( stripID<1) {
;
} else if (stripID<20) {
xa=1; ya=2.;
xb=20; yb=4.;
} else if (stripID<250) {
xa=20; ya=4.;
xb=250; yb=6.;
} else {
xa=250; ya=6.;
xb=290; yb=9.;
}
y=ya+(yb-ya)/(xb-xa) *(stripID -xa);
Float_t gain = 8.0*y;
// ----<<<<< temporary kludge for gains >>>>>----
//
(m_EEmcDbase -> byIndex+index) -> gain = gain;
}
}
}
}
/////////////////////////////////////////////////////////////////////////////
void EEsmdCalibration::newDbRatio ( Int_t ifirst, Int_t ilast )
{
// Calculates new database values based on the histograms
// of the event-by-event ratios of data yield / fit.
assert(m_EEmcDbase);
TH1F *nentr;
TH1F *hratio;
TH1F *hwidth;
Int_t first = 1;
if ( first ) {
cd();
TDirectory *sub = new TDirectory("newDbRatio","New gains from ratios of data yields to fits");
sub -> cd();
first = 0;
nentr = new TH1F("nentr","Number of entries in ratio histograms",100,0.,400.);
hratio = new TH1F("hratio","Distrib of yield-to-fit ratios",100,0.,2.);
hwidth = new TH1F("hwidthratio","Width of dist. of yield-to-fit ratios",100,0.,2.);
sub -> cd("..");
cd("..");
}
//-- Sector loop --
for ( Int_t i = ifirst; i <= ilast; i++ ) {
for ( Int_t j = 0; j < 2; j++ ) {
Float_t old_sum = 0.;
Float_t new_sum = 0.;
for ( Int_t k = 0; k < 288; k++ ) {
// Name of the strip
TString strip;
if ( i+1 < 10 ) strip += "0"; strip += i+1;
strip += (j) ? "V" : "U";
if ( k+1<100 ) strip += "0"; if ( k+1<10 ) strip += "0"; strip += k+1;
// Index of the database entry
Int_t index = EEname2Index( strip );
// Database entry
const EEmcDbItem *x = m_EEmcDbase -> getByIndex(index);
// Verify that the DB entry is Ok
if ( !x->fail && !x->stat )
old_sum += x -> gain;
else
continue;
// Do not even attempt to fit if too few entries
Int_t nentries = (Int_t)m_YieldToFit[i][j][k] -> GetEntries();
nentr -> Fill( (Float_t)nentries );
if ( nentries < 30 ) {
std::cout << "Abort correction, insufficient entries, " << x -> name << std::endl;
continue;
}
std::cout << std::endl
<< "############################################"
<< std::endl
<< x->name
<< std::endl;
// Fit the ratio dist. for this strip to a gaussian
TF1 fit("fit","gaus",0.3,1.5);
m_YieldToFit[i][j][k] -> Fit(&fit,"R");
Float_t ratio = fit.GetParameter(1);
Float_t width = fit.GetParameter(2);
Float_t eratio = fit.GetParError(1);
Float_t ewidth = fit.GetParError(2);
//--
//-- Apply some sanity checks before we correct gains
//--
//--
//-- Do not make too large of a correction
//--
hratio -> Fill( ratio );
if ( ratio < 0.25 || ratio > 4.0 ) {
std::cout << "Abort correction, too large, " << x -> name << std::endl;
continue;
}
//--
//-- If the width is too wide (need to define "too wide")
//-- skip it.
//--
hwidth -> Fill( width );
if ( width > 0.5 ) {
continue;
}
//--
//-- Increase or decrease the ratio by 5% for each iteration
//-- to "slow down" the convergence
//--
if ( ratio < 1. )
ratio *= 1.05;
else if ( ratio > 1. )
ratio /= 1.05;
//--
//-- Now correct the database
//--
(m_EEmcDbase -> byIndex+index) -> gain *= ratio;
new_sum += m_EEmcDbase -> getByIndex(index) -> gain;
}//-- Loop over strips
for ( Int_t k = 0; k < 288; k++ ) {
TString strip;
if ( i+1 < 10 ) strip += "0"; strip += i+1;
strip += (j) ? "V" : "U";
if ( k+1<100 ) strip += "0"; if ( k+1<10 ) strip += "0"; strip += k+1;
// Index of the database entry
Int_t index = EEname2Index( strip );
(m_EEmcDbase -> byIndex+index) -> gain *= (old_sum / new_sum);
}//-- Loop over strips (renormalization)
}
}
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//t_t EEsmdCalibration::isolationCut ( EEezCluster *cluster )
Int_t EEsmdCalibration::isolationCut ( Int_t icluster )
{
// Determine if the given cluster is isolated in the sector,
// in the sense that it is unlikely that neighboring towers
// contain hits which may contaminate the SMD strips in the
// seed tower.
//
// A return value of 1 indicates the cluster FAILED the
// isolation cut.
//
// This will be a simple isolation cut. Any tower within +/- 3
// etabins which is occupied will trigger the cut (any tower
// not within the cluster that is.)
//
EEezCluster *cluster = m_EEezClust -> getClusterPtr(icluster);
static Int_t first=1;
static TH1F *h_iso_energy;
static TH1F *h_iso_towers;
static TH1F *h_iso_emaxim;
static TH2F *h_iso_xy;
static TH1F *h_iso_dndeta;
static TH1F *h_iso_strips[2];
static TH1F *h_iso_nmips[2];
static TH1F *h_iso_ecluster;
static TH2F *h_iso_frac9_vs_eiso[2];
static TH1F *h_pass_iso_energy;
static TH1F *h_pass_iso_towers;
static TH1F *h_pass_iso_emaxim;
static TH2F *h_pass_iso_xy;
static TH1F *h_pass_iso_dndeta;
static TH1F *h_pass_iso_strips[2];
static TH1F *h_pass_iso_nmips[2];
static TH1F *h_pass_iso_ecluster;
static TH2F *h_pass_iso_frac9_vs_eiso[2];
if ( first ) {
cd();
TDirectory *sub = new TDirectory("isolationCut","Histogams showing various quantities before isolation cuts");
sub -> cd();
first = 0;
h_iso_energy = new TH1F("h_iso_energy", "Energy sum of towers in isolation cut of cluster",50,0.,2.);
h_iso_towers = new TH1F("h_iso_towers", "Number of towers within isolation cut of cluster",20,0.,20.);
h_iso_emaxim = new TH1F("h_iso_emaxim", "Maximum energy of tower within isolation cut of cluster",50,0.,2.);
h_iso_xy = new TH2F("h_iso_xy", "Tower-weighted Y vs X of clusters",100,-175.,175.,100,-175.,175.);
h_iso_dndeta = new TH1F("h_iso_dndeta", "Number of clusters vs eta-bin",12,0.,12.);
h_iso_strips[0] = new TH1F("h_iso_stripsU","Number of hit strips beneath seed tower versus strip index U planes",288,0.,288.);
h_iso_strips[1] = new TH1F("h_iso_stripsV","Number of hit strips beneath seed tower versus strip index V planes",288,0.,288.);
h_iso_nmips[0] = new TH1F("h_iso_nmipsU","Number of SMD mips beneath seed tower versus strip index U planes",288,0.,288.);
h_iso_nmips[1] = new TH1F("h_iso_nmipsV","Number of SMD mips beneath seed tower versus strip index V planes",288,0.,288.);
h_iso_ecluster = new TH1F("h_iso_ecluster","Energy of all cluster",50,0.,10.);
h_iso_frac9_vs_eiso[0] = new TH2F("h_iso_frac9_vs_eisoU","Fractional number of mips in 9 strips adjacent to max strip versus E in iso zone",100,0.,0.8,50,0.,1.);
h_iso_frac9_vs_eiso[1] = new TH2F("h_iso_frac9_vs_eisoV","Fractional number of mips in 9 strips adjacent to max strip versus E in iso zone",100,0.,0.8,50,0.,1.);
/* sub -> cd("..");
sub = new TDirectory("passIsolationCut","Histograms showing various quantities after isolation cuts");
sub -> cd();*/
h_pass_iso_energy = new TH1F("h_pass_iso_energy", "Energy sum of towers in isolation cut of cluster",50,0.,2.);
h_pass_iso_towers = new TH1F("h_pass_iso_towers", "Number of towers within isolation cut of cluster",20,0.,20.);
h_pass_iso_emaxim = new TH1F("h_pass_iso_emaxim", "Maximum energy of tower within isolation cut of cluster",50,0.,2.);
h_pass_iso_xy = new TH2F("h_pass_iso_xy", "Tower-weighted Y vs X of clusters",100,-175.,175.,100,-175.,175.);
h_pass_iso_dndeta = new TH1F("h_pass_iso_dndeta", "Number of clusters vs eta-bin",12,0.,12.);
h_pass_iso_strips[0] = new TH1F("h_pass_iso_stripsU","Number of hit strips beneath seed tower versus strip index U planes",288,0.,288.);
h_pass_iso_strips[1] = new TH1F("h_pass_iso_stripsV","Number of hit strips beneath seed tower versus strip index V planes",288,0.,288.);
h_pass_iso_nmips[0] = new TH1F("h_pass_iso_nmipsU","Number of SMD mips beneath seed tower versus strip index U planes",288,0.,288.);
h_pass_iso_nmips[1] = new TH1F("h_pass_iso_nmipsV","Number of SMD mips beneath seed tower versus strip index V planes",288,0.,288.);
h_pass_iso_ecluster = new TH1F("h_pass_iso_ecluster","Energy of all cluster",50,0.,10.);
h_pass_iso_frac9_vs_eiso[0] = new TH2F("h_pass_iso_frac9_vs_eisoU","Fractional number of mips in 9 strips adjacent to max strip versus E in iso zone",100,0.,0.8,50,0.,1.);
h_pass_iso_frac9_vs_eiso[1] = new TH2F("h_pass_iso_frac9_vs_eisoV","Fractional number of mips in 9 strips adjacent to max strip versus E in iso zone",100,0.,0.8,50,0.,1.);
sub -> cd("..");
cd("..");
}
//--
//-- Cluster energy
//--
h_iso_ecluster -> Fill ( cluster -> getEnergy() );
//--
//-- x vs y of clusters (tower-weighted)
//--
TVector3 momentum = cluster -> getMomentum();
TVector3 direction = momentum.Unit();
h_iso_xy->Fill( (kEEmcZSMD*direction).X(), (kEEmcZSMD*direction).Y() );
//--
//-- Eta bin of the cluster
//--
Int_t sec = cluster -> getSeedTower() -> getSector();
Int_t sub = cluster -> getSeedTower() -> getSubSector();
Int_t eta = cluster -> getSeedTower() -> getEtabin();
h_iso_dndeta -> Fill ( (Float_t)eta + 0.5 );
//--
//-- Number of his strips beneath the seed tower
//--
EEezTower *seed = cluster -> getSeedTower();
for ( Int_t i = 0; i < seed -> getNUStrips(); i++ ) {
EEezStrip *strip = seed -> getUStrip(i);
if ( strip -> getEnergy() > 0. ) {
h_iso_strips[0] -> Fill( strip -> getIndex() );
h_iso_nmips[0] -> Fill( strip -> getIndex()+0.5, strip -> getNumMips() );
}
}
for ( Int_t i = 0; i < seed -> getNVStrips(); i++ ) {
EEezStrip *strip = seed -> getVStrip(i);
if ( strip -> getEnergy() > 0. ) {
h_iso_strips[1] -> Fill( strip -> getIndex() );
h_iso_nmips[1] -> Fill( strip -> getIndex()+0.5, strip -> getNumMips() );
}
}
Int_t min = TMath::Max(0, eta-3);
Int_t max = TMath::Min(12,eta+3);
Float_t esum = 0.;
Int_t stat = 0;
Float_t emax = 0.;
Int_t ntow = 0; // Number of struck towers in isolation zone
//--
//-- Kills the event (stat != 0) if any one tower has > 0.5 GeV energy
//-- deposit within the +/- 3 bins in eta of the tower
//--
for ( Int_t myeta = min; myeta <= max; myeta++ ) {
for ( Int_t mysub = 0; mysub < 5; mysub++ ) {
Int_t deta = TMath::Abs(myeta-eta);
Int_t dsub = TMath::Abs(mysub-sub);
// Ignore towers in cluster and towers at the same
// eta and phi
if ( deta <= 1 && dsub <= 1 ) continue;
EEezTower *tower = m_EEezAnaly -> getTower(sec,mysub,myeta);
Float_t etow = tower -> getEnergy();
if ( etow <= 0. ) continue;
if ( etow > 0.5 ) stat++;
ntow++;
esum += etow;
if ( etow > emax ) emax = etow;
h_iso_energy -> Fill( etow );
}
}
//--
//-- Fractional energy deposited near SMD's max strip vs iso region
//--
TH1F *profileU = m_EEsmdProf -> getProfileU(icluster);
TH1F *profileV = m_EEsmdProf -> getProfileV(icluster);
Int_t maxU = profileU -> GetMaximumBin();
Int_t maxV = profileV -> GetMaximumBin();
Float_t sumU = profileU -> Integral();
Float_t sumV = profileV -> Integral();
Float_t sumU9 = profileU -> Integral(maxU-4,maxU+4);
Float_t sumV9 = profileV -> Integral(maxV-4,maxV+4);
if ( sumU > 0. ) {
h_iso_frac9_vs_eiso[0] -> Fill(esum, sumU9/sumU);
if ( stat == 0 )
h_pass_iso_frac9_vs_eiso[0] -> Fill(esum, sumU9/sumU);
}
if ( sumV > 0. ) {
h_iso_frac9_vs_eiso[1] -> Fill(esum, sumV9/sumV);
if ( stat == 0 )
h_pass_iso_frac9_vs_eiso[1] -> Fill(esum, sumV9/sumV);
}
//--
//--
//--
if ( emax > 0. )
h_iso_emaxim -> Fill( emax );
h_iso_towers -> Fill( (Float_t)ntow + 0.5 );
//--
//-- Fill histograms for events which pass isolation cuts
//--
if ( stat == 0 ) {
//$$$std::cout << "Isolation cut passed" << std::endl;
h_pass_iso_ecluster -> Fill( cluster -> getEnergy() );
h_pass_iso_xy->Fill( (kEEmcZSMD*direction).X(), (kEEmcZSMD*direction).Y() );
h_pass_iso_dndeta -> Fill ( (Float_t)eta + 0.5 );
//--
//-- Number of his strips beneath the seed tower
//--
EEezTower *seed = cluster -> getSeedTower();
for ( Int_t i = 0; i < seed -> getNUStrips(); i++ ) {
EEezStrip *strip = seed -> getUStrip(i);
if ( strip -> getEnergy() > 0. ) {
h_pass_iso_strips[0] -> Fill( strip -> getIndex() );
h_pass_iso_nmips[0] -> Fill( strip -> getIndex()+0.5, strip -> getNumMips() );
}
}
for ( Int_t i = 0; i < seed -> getNVStrips(); i++ ) {
EEezStrip *strip = seed -> getVStrip(i);
if ( strip -> getEnergy() > 0. ) {
h_pass_iso_strips[1] -> Fill( strip -> getIndex() );
h_pass_iso_nmips[1] -> Fill( strip -> getIndex()+0.5, strip -> getNumMips() );
}
}
for ( Int_t myeta = min; myeta <= max; myeta++ ) {
for ( Int_t mysub = 0; mysub < 5; mysub++ ) {
Int_t deta = TMath::Abs(myeta-eta);
Int_t dsub = TMath::Abs(mysub-sub);
// Ignore towers in cluster and towers at the same
// eta and phi
if ( deta <= 1 && dsub <= 1 && ( deta!=0 || dsub!=0 ) ) continue;
EEezTower *tower = m_EEezAnaly -> getTower(sec,mysub,myeta);
Float_t etow = tower -> getEnergy();
if ( etow <= 0. ) continue;
h_pass_iso_energy -> Fill( etow );
}
}
if ( emax > 0. )
h_pass_iso_emaxim -> Fill( emax );
h_pass_iso_towers -> Fill( (Float_t)ntow + 0.5 );
}
return stat;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEsmdCalibration::scaleGains( Float_t scale, Int_t ifirst, Int_t ilast )
{
// Apply a one-time normalization to the gains for all planes
for ( Int_t i = ifirst; i <= ilast; i++ ) {
TString sec = "";
if ( (i+1)<10 ) sec+= "0";
sec += (i+1);
for ( Int_t j = 0; j < 2; j++ ) {
TString plane = sec;
plane += ( j ) ? "V" : "U";
TString histo = "gshift"; histo += plane;
cd();
TH1F *h = (TH1F *) Get(histo);
cd("..");
std::cout << "Processing plane = " << plane << std::endl << std::flush;
Float_t oldSum = 0.;
Float_t newSum = 0.;
for ( Int_t k = 0; k < 288; k++ ) {
TString name = plane;
if ( (k+1)<100 ) name += "0";
if ( (k+1)<10 ) name += "0";
name += (k+1);
std::cout << name << " " << std::flush;
Int_t index = EEname2Index( name );
std::cout << index << " " << std::flush;
const EEmcDbItem *x = m_EEmcDbase -> getByIndex(index);
std::cout << x -> gain << " ";
Float_t gain = scale * x -> gain;
(m_EEmcDbase -> byIndex+index) -> gain = gain;
std::cout << x -> gain << std::endl;
}
}
}
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
Int_t EEsmdCalibration::fitQAcut ( Int_t icluster, Int_t iuv )
{
// Perform simple QA cuts on the profile histogram(s). A return
// value of "1" indicates the cut was triggered and the event
// should be excluded from further analysis.
// Book histograms for QA studies (inefficient, but make the code
// slightly more readable... good considering how bloated this code
// is becoming...).
static TH1F *h_fitqa_errmean[2];
static TH1F *h_fitqa_errwidth[2];
static TH1F *h_fitqa_erryield[2];
static TH1F *h_fitqa_yield[2];
static TH1F *h_fitqa_width[2];
static TH2F *h_fitqa_yield_vs_ec[2]; //----<<<< need to fill, and create pass
static TH1F *h_fitqa_chi2[2];
static TH2F *h_fitqa_chi2_vs_strip[2];
static TH2F *h_pass_fitqa_xy[2];
static TH1F *h_pass_fitqa_dndeta[2];
static TH1F *h_pass_fitqa_strips[2];
static TH1F *h_pass_fitqa_nmips[2];
static TH1F *h_pass_fitqa_ecluster[2];
static TH1F *h_pass_fitqa_chi2[2];
static TH2F *h_pass_fitqa_chi2_vs_strip[2];
TDirectory *d_fail_errmean;
Int_t i_fail_errmean = 0;
TDirectory *d_fail_width;
Int_t i_fail_width = 0;
TDirectory *d_pass_fitqa;
Int_t i_pass_fitqa = 0;
static Int_t first = 1;
if ( first ) {
first = 0;
cd();
TDirectory *sub = new TDirectory("fitQAcut","Histograms of quantities which pass QA cuts on fits");
sub -> cd();
d_fail_errmean = new TDirectory("fail_errmean","SMD profile histograms which fail QA cut on mean uncertainty");
d_fail_width = new TDirectory("fail_width", "SMD profile histograms which fail QA cut on width");
d_pass_fitqa = new TDirectory("pass_fitqa", "SMD profile histograms which pass QA cuts");
//-- Things before the fits
h_fitqa_erryield[0]= new TH1F("h_fitqa_erryieldU", "Error on fit yield (fractional)", 100,0.,1. );
h_fitqa_errmean[0] = new TH1F("h_fitqa_errmeanU", "Error on fit mean", 100, 0., 2.);
h_fitqa_errwidth[0]= new TH1F("h_fitqa_errwidthU", "Error on fit width", 100, 0., 2.);
h_fitqa_yield[0] = new TH1F("h_fitqa_yieldU", "Yield of mips",100,0.,200.);
h_fitqa_width[0] = new TH1F("h_fitqa_widthU", "Width", 50, 0., 5. );
h_fitqa_chi2[0] = new TH1F("h_fitqa_chi2U", "#chi^2 / ndf distribution",100,0.,25.);
h_fitqa_chi2_vs_strip[0] = new TH2F("h_fitqa_chi2_vs_stripU","#chi^2 vs strip number",288,0.,288.,100,0.,25.);
h_fitqa_erryield[1]= new TH1F("h_fitqa_erryieldV", "Error on fit yield (fractional)", 100,0.,1. );
h_fitqa_errmean[1] = new TH1F("h_fitqa_errmeanV", "Error on fit mean", 100, 0., 2.);
h_fitqa_errwidth[1]= new TH1F("h_fitqa_errwidthV", "Error on fit width", 100, 0., 2.);
h_fitqa_yield[1] = new TH1F("h_fitqa_yieldV", "Yield of mips",100,0.,200.);
h_fitqa_width[1] = new TH1F("h_fitqa_widthV", "Width", 50, 0., 5. );
h_fitqa_chi2[1] = new TH1F("h_fitqa_chi2V", "#chi^2 / ndf distribution",100,0.,25.);
h_fitqa_chi2_vs_strip[1] = new TH2F("h_fitqa_chi2_vs_stripV","#chi^2 vs strip number",288,0.,288.,100,0.,25.);
//-- Quantities which we like to monitor as cuts are applied
h_pass_fitqa_xy[0] = new TH2F("h_pass_fitqa_xyU", "Tower-weighted Y vs X of clusters",100,-175.,175.,100,-175.,175.);
h_pass_fitqa_dndeta[0] = new TH1F("h_pass_fitqa_dndetaU", "Number of clusters vs eta-bin",12,0.,12.);
h_pass_fitqa_xy[1] = new TH2F("h_pass_fitqa_xyV", "Tower-weighted Y vs X of clusters",100,-175.,175.,100,-175.,175.);
h_pass_fitqa_dndeta[1] = new TH1F("h_pass_fitqa_dndetaV", "Number of clusters vs eta-bin",12,0.,12.);
h_pass_fitqa_strips[0] = new TH1F("h_pass_fitqa_stripsU","Number of hit strips beneath seed tower versus strip index U planes",288,0.,288.);
h_pass_fitqa_strips[1] = new TH1F("h_pass_fitqa_stripsV","Number of hit strips beneath seed tower versus strip index V planes",288,0.,288.);
h_pass_fitqa_nmips[0] = new TH1F("h_pass_fitqa_nmipsU","Number of SMD mips beneath seed tower versus strip index U planes",288,0.,288.);
h_pass_fitqa_nmips[1] = new TH1F("h_pass_fitqa_nmipsV","Number of SMD mips beneath seed tower versus strip index V planes",288,0.,288.);
h_pass_fitqa_ecluster[0] = new TH1F("h_pass_fitqa_eclusterU","Energy of all cluster",50,0.,10.);
h_pass_fitqa_ecluster[1] = new TH1F("h_pass_fitqa_eclusterV","Energy of all cluster",50,0.,10.);
//-- Filled for U, even if V fails (and vice versa?)
h_pass_fitqa_chi2[0] = new TH1F("h_pass_fitqa_chi2U", "#chi^2 / ndf distribution",100,0.,25.);
h_pass_fitqa_chi2_vs_strip[0] = new TH2F("h_pass_fitqa_chi2_vs_stripU","#chi^2 vs strip number",288,0.,288.,100,0.,25.);
h_pass_fitqa_chi2[1] = new TH1F("h_pass_fitqa_chi2V", "#chi^2 / ndf distribution",100,0.,25.);
h_pass_fitqa_chi2_vs_strip[1] = new TH2F("h_pass_fitqa_chi2_vs_stripV","#chi^2 vs strip number",288,0.,288.,100,0.,25.);
sub -> cd("..");
cd("..");
}
//-- Get the appropriate profile histogram
TH1F *profile = ( iuv ) ?
m_EEsmdProf -> getProfileV(icluster) :
m_EEsmdProf -> getProfileU(icluster) ;
//-- Get the corresponding cluster
EEezCluster *cluster = m_EEezClust -> getClusterPtr(icluster);
//-- And a pointer to its fit
TF1 *fit = (TF1 *)profile -> GetListOfFunctions() -> At(0);
Float_t yield = fit -> GetParameter(0);
Float_t mean = fit -> GetParameter(1);
Float_t width = fit -> GetParameter(2);
Float_t eyield = fit -> GetParError(0);
Float_t emean = fit -> GetParError(1);
Float_t ewidth = fit -> GetParError(2);
Float_t chi2 = fit -> GetChisquare();
Float_t ndf = fit -> GetNDF();
Int_t mean_bin = profile -> FindBin( mean );
Int_t count = 0;
Float_t sumstrips5 = 0.;
Float_t sumstrips9 = 0.;
//-- Cluster energy
Float_t energy = cluster -> getEnergy();
//-- Chi-squared distributions
if ( ndf > 0. ) {
h_fitqa_chi2[iuv] -> Fill ( chi2 / ndf );
h_fitqa_chi2_vs_strip[iuv] -> Fill ( profile -> FindBin(mean) - 0.5, chi2/ndf );
}
//--
//-- Apply sanity checks on the parameters of the fits.
//--
if ( yield > 0. )
h_fitqa_erryield[iuv] -> Fill( eyield / yield );
h_fitqa_errmean[iuv] -> Fill( emean );
h_fitqa_errwidth[iuv] -> Fill( ewidth );
h_fitqa_yield[iuv] -> Fill ( yield );
if ( emean > 0.2 ) { // return 1; // err mean > 1 strip, punt
// if ( i_fail_errmean < 50 ) {
// TString hname = profile -> GetName();
// hname += "_fail_errmean_";
// hname += i_fail_errmean++;
// d_fail_errmean -> cd();
// profile -> Clone(hname);
// d_fail_errmean -> cd("..");
// }
return 1;
}
h_fitqa_width[iuv] -> Fill ( width );
if ( width < 0.5 || width > 2.0 ) {
// if ( i_fail_width < 50 ) {
// TString hname = profile -> GetName();
// hname += "_fail_width_";
// hname += i_fail_width++;
// d_fail_width -> cd();
// profile -> Clone(hname);
// d_fail_width -> cd("..");
// }
return 1;
}
//-- Chi-squared distributions
if ( ndf > 0. ) {
h_pass_fitqa_chi2[iuv] -> Fill ( chi2 / ndf );
h_pass_fitqa_chi2_vs_strip[iuv] -> Fill ( profile -> FindBin(mean) - 0.5, chi2/ndf );
}
//--
//-- Fill passed event quantities
//--
TVector3 momentum = cluster -> getMomentum();
TVector3 direction = momentum.Unit();
h_pass_fitqa_xy[iuv] ->Fill( (kEEmcZSMD*direction).X(), (kEEmcZSMD*direction).Y() );
Int_t eta = cluster -> getSeedTower() -> getEtabin();
h_pass_fitqa_dndeta[iuv] -> Fill ( (Float_t)eta + 0.5 );
//--
//-- Number of his strips beneath the seed tower
//--
EEezTower *seed = cluster -> getSeedTower();
if ( iuv == 0 )
for ( Int_t i = 0; i < seed -> getNUStrips(); i++ ) {
EEezStrip *strip = seed -> getUStrip(i);
if ( strip -> getEnergy() > 0. ) {
h_pass_fitqa_strips[0] -> Fill( strip -> getIndex() );
h_pass_fitqa_nmips[0] -> Fill( strip -> getIndex()+0.5, strip -> getNumMips() );
}
}
else
for ( Int_t i = 0; i < seed -> getNVStrips(); i++ ) {
EEezStrip *strip = seed -> getVStrip(i);
if ( strip -> getEnergy() > 0. ) {
h_pass_fitqa_strips[1] -> Fill( strip -> getIndex() );
h_pass_fitqa_nmips[1] -> Fill( strip -> getIndex()+0.5, strip -> getNumMips() );
}
}
h_pass_fitqa_ecluster[iuv] -> Fill ( cluster -> getEnergy() );
//--
//-- Save first 50 SMD profiles which pass QA cuts on fits
//--
if ( i_pass_fitqa < 50 ) {
// TString hname = profile -> GetName();
// hname += "_pass_fitqa_";
// hname += i_pass_fitqa++;
// d_pass_fitqa -> cd();
// profile -> Clone(hname);
// d_pass_fitqa -> cd("..");
}
return 0;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
Int_t EEsmdCalibration::profileQAcut( Int_t icluster, Int_t iuv )
{
// Perform simple QA analysis of the SMD distributions beneath
// an identified, isolated shower
assert( iuv==0 || iuv==1 );
static TH1F *h_profqa_nof5[2];
static TH1F *h_profqa_frac5to9[2];
static TH1F *h_profqa_frac9toAll[2];
static TH2F *h_pass_profqa_xy[2];
static TH1F *h_pass_profqa_dndeta[2];
static TH1F *h_pass_profqa_strips[2];
static TH1F *h_pass_profqa_nmips[2];
static TH1F *h_pass_profqa_ecluster[2];
static TDirectory *d_fail_count;
static Int_t i_fail_count = 0;
static TDirectory *d_fail_5to9;
static Int_t i_fail_5to9 = 0;
static TDirectory *d_pass_profqa;
static Int_t i_pass_profqa = 0;
static Int_t first = 1;
if ( first ) {
cd();
TDirectory *sub = new TDirectory("profileQAcut","Histograms of quantities which pass rudimentary strip-wise QA cuts");
sub -> cd();
d_fail_count = new TDirectory("fail_count","Profile histograms which fail the 4/5 strips criteria");
d_fail_5to9 = new TDirectory("fail_5to9", "Profile histograms which fail the 5-strip-sum to 9-strip-sum");
d_pass_profqa = new TDirectory("pass_profqa", "Profile histograms which pass rudimentary profile QA");
first = 0;
//-- Quantities we cut on
h_profqa_nof5[0] = new TH1F("h_profqa_nof5U", "Number of strips within #pm 2 strips of maximum bin ", 10, 0., 10. );
h_profqa_frac5to9[0] = new TH1F("h_profqa_frac5to9U","Ratio of number of mips deposited w/in +/- 2 strip of max bin to that deposited w/in +/- 4 strips", 50, 0., 1.);
h_profqa_frac9toAll[0] = new TH1F("h_profqa_frac9toAllU","Ratio of number of mips deposited w/in +/- 2 strips of max bin to that deposited w/in full range under seed tower (no cut yet)", 50, 0., 1.);
h_profqa_nof5[1] = new TH1F("h_profqa_nof5V", "Number of strips within #pm 2 strips of maximum bin ", 10, 0., 10. );
h_profqa_frac5to9[1] = new TH1F("h_profqa_frac5to9V","Ratio of number of mips deposited w/in +/- 2 strip of max bin to that deposited w/in +/- 4 strips", 50, 0., 1.);
h_profqa_frac9toAll[1] = new TH1F("h_profqa_frac9toAllV","Ratio of number of mips deposited w/in +/- 2 strips of max bin to that deposited w/in full range under seed tower (no cut yet)", 50, 0., 1.);
//-- Quantities we like to monitor from cut-to-cut... in this case, things which passed
h_pass_profqa_xy[0] = new TH2F("h_pass_profqa_xyU", "Tower-weighted Y vs X of clusters",100,-175.,175.,100,-175.,175.);
h_pass_profqa_dndeta[0] = new TH1F("h_pass_profqa_dndetaU", "Number of clusters vs eta-bin",12,0.,12.);
h_pass_profqa_xy[1] = new TH2F("h_pass_profqa_xyV", "Tower-weighted Y vs X of clusters",100,-175.,175.,100,-175.,175.);
h_pass_profqa_dndeta[1] = new TH1F("h_pass_profqa_dndetaV", "Number of clusters vs eta-bin",12,0.,12.);
h_pass_profqa_strips[0] = new TH1F("h_pass_profqa_stripsU","Number of hit strips beneath seed tower versus strip index U planes",288,0.,288.);
h_pass_profqa_strips[1] = new TH1F("h_pass_profqa_stripsV","Number of hit strips beneath seed tower versus strip index V planes",288,0.,288.);
h_pass_profqa_nmips[0] = new TH1F("h_pass_profqa_nmipsU","Number of SMD mips beneath seed tower versus strip index U planes",288,0.,288.);
h_pass_profqa_nmips[1] = new TH1F("h_pass_profqa_nmipsV","Number of SMD mips beneath seed tower versus strip index V planes",288,0.,288.);
h_pass_profqa_ecluster[0] = new TH1F("h_pass_profqa_eclusterU","Energy of all cluster",50,0.,10.);
h_pass_profqa_ecluster [1] = new TH1F("h_pass_profqa_eclusterV","Energy of all cluster",50,0.,10.);
sub -> cd("..");
cd("..");
}
//-- Get the appropriate profile histogram
TH1F *profile = ( iuv ) ?
m_EEsmdProf -> getProfileV(icluster) :
m_EEsmdProf -> getProfileU(icluster) ;
EEezCluster *cluster =
m_EEezClust -> getClusterPtr(icluster);
//-- The fitting algorithm currently attempts to fit the histogram
//-- by guessing that the mean is near the strip with the largest
//-- energy response.
//-- Our first QA test requires that of the five consecutive strips
//-- including the maximum strip, we have four of the five strips
//-- showing a nonzero energy response.
Int_t mean_bin = profile -> GetMaximumBin();
Int_t count = 0;
Float_t sumstrips5 = 0.;
Float_t sumstrips9 = 0.;
Float_t sumstripsAll = profile -> Integral();
Int_t stat = 0;
for ( Int_t i = mean_bin - 2; i <= mean_bin +2; i++ ) {
Float_t nm=profile -> GetBinContent( i );
if ( nm > 0.0 ) {
count++;
sumstrips5 += nm;
}
}
h_profqa_nof5[iuv] -> Fill ( count );
if ( count < 4 ) {
TString hname = profile -> GetName();
hname += "_failed_count_";
hname += i_fail_count++;
if ( i_fail_count < 50 ) {
d_fail_count -> cd();
profile -> Clone(hname);
d_fail_count -> cd("..");
}
return 1;
}
//--
//-- Next, sum up number of mips in adjacent 4 strips on either side
//--
for ( Int_t i = mean_bin - 4; i <= mean_bin + 4; i++ ) {
sumstrips9 += profile -> GetBinContent( i );
}
if ( sumstrips9 <= 0. ) return 1;
h_profqa_frac5to9[iuv] -> Fill ( sumstrips5 / sumstrips9 );
if ( sumstrips5 / sumstrips9 < 0.6 ) {
TString hname = profile -> GetName();
hname += "_failed_5to9_";
hname += i_fail_5to9;
if ( i_fail_5to9 < 50 ) {
d_fail_5to9 -> cd();
profile -> Clone(hname);
d_fail_5to9 -> cd("..");
}
return 1;
}
if ( sumstripsAll > 0. ) h_profqa_frac9toAll[iuv] -> Fill ( sumstrips9 / sumstripsAll );
//--
//-- Fill histograms for "passed" quantities.
//--
TVector3 momentum = cluster -> getMomentum();
TVector3 direction = momentum.Unit();
h_pass_profqa_xy[iuv] ->Fill( (kEEmcZSMD*direction).X(), (kEEmcZSMD*direction).Y() );
Int_t eta = cluster -> getSeedTower() -> getEtabin();
h_pass_profqa_dndeta[iuv] -> Fill ( (Float_t)eta + 0.5 );
//--
//-- Number of his strips beneath the seed tower
//--
EEezTower *seed = cluster -> getSeedTower();
if ( iuv == 0 )
for ( Int_t i = 0; i < seed -> getNUStrips(); i++ ) {
EEezStrip *strip = seed -> getUStrip(i);
if ( strip -> getEnergy() > 0. ) {
h_pass_profqa_strips[0] -> Fill( strip -> getIndex() );
h_pass_profqa_nmips[0] -> Fill( strip -> getIndex()+0.5, strip -> getNumMips() );
}
}
else
for ( Int_t i = 0; i < seed -> getNVStrips(); i++ ) {
EEezStrip *strip = seed -> getVStrip(i);
if ( strip -> getEnergy() > 0. ) {
h_pass_profqa_strips[1] -> Fill( strip -> getIndex() );
h_pass_profqa_nmips[1] -> Fill( strip -> getIndex()+0.5, strip -> getNumMips() );
}
}
h_pass_profqa_ecluster[iuv] -> Fill ( cluster -> getEnergy() );
//--
//-- Save this profile histogram
//--
if ( i_pass_profqa < 50 ) {
d_pass_profqa -> cd();
TString hname = profile -> GetName();
hname += "_passed_profqa_";
hname += i_pass_profqa++;
profile -> Clone(hname);
d_pass_profqa -> cd("..");
}
return 0;
}
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.