#include "EEezPi0Analysis.h"
#include <TString.h>
#include <TMath.h>
#include "StEEmcUtil/EEfeeRaw/EEstarTrig.h"
#include "EEezAnalysis.h"
#include <iostream>
#include <algorithm> // provides random_shuffle
/////////////////////////////////////////////////////////////////////////////
ClassImp(EEezPi0Analysis);
EEezPi0Analysis::EEezPi0Analysis( const Char_t *name, Int_t trigger_id )
{
// The constructor. Takes a name (for the directory) and
// trigger ID. Events with different trigger ID's are not
// processed, since they would have potentially very different
// topologies.
m_TriggerId = trigger_id;
m_Name = name;
m_Directory = new TDirectory( m_Name, "EEmc pi0 analysis" );
setCtbThreshold();
m_Random = new TRandom(); // with default seed... no method to set seed,
// as we want reproducability between any two
// runs (random numbers used to mix backgrounds
// only.
m_EEmcGeom = new EEmcGeomSimple();
m_CtbMinMultiplicity = 0.;
m_CtbMaxMultiplicity = 240.;
m_CtbMinAdcSum = 0;
m_CtbMaxAdcSum = 240 * 256;
m_MinMass = 0.135 - 3.0*0.03;
m_MaxMass = 0.135 + 3.0*0.03;
}
/////////////////////////////////////////////////////////////////////////////
Int_t EEezPi0Analysis::Init()
{
// Initialize pi0 finder.
InitHistograms();
return 1;
}
void EEezPi0Analysis::InitHistograms()
{
// Book an ever-increasing multitude of histograms
m_Directory -> cd();
TList *mylist = new TList(); // temp list;
//
// Initialize histograms for kinematic quantities
//
m_Mass = new TH1F("m_Mass", "Two-cluster invariant mass distribution", 100,0.,1. );
mylist -> Add(m_Mass);
m_Mass1 = new TH1F("m_Mass1", "Two-cluster invariant mass distribution", 200,0.,1. );
mylist -> Add(m_Mass1);
m_Mass2 = new TH1F("m_Mass2", "Two-cluster invariant mass distribution", 300,0.,1. );
mylist -> Add(m_Mass2);
m_Energy = new TH1F("m_Energy", "Two-cluster energy", 100, 0., 20.);
mylist -> Add(m_Energy);
m_Z = new TH1F("m_Z", "Two-cluster Z = |E#_1 - E#_2| / (E#_1+E#_2)", 100, 0.,1. );
mylist -> Add(m_Z);
m_Eta = new TH1F("m_Eta", "Two-cluster #eta distribution", 100, 1., 2.);
mylist -> Add(m_Eta);
m_Phi = new TH1F("m_Phi", "Two-cluster #phi distribution", 100, -TMath::Pi(), TMath::Pi() );
mylist -> Add(m_Phi);
// Trigger
m_CtbAdc = new TH1F("m_CtbAdc", "CTB adc response, for 240 individual towers", 256,0.,256.);
//
// Kinematic quantities, binned by etabin, phibin and tower
//
for ( Int_t ieta = 0; ieta < 12; ieta++ ) {
TString index = "["; index += ieta; index += "]";
TString bin = " #eta bin number "; bin += (ieta+1);
m_Mass_DEta[ieta] = new TH1F( TString("m_Mass_DEta") +index, TString("Two-cluster invariant mass")+bin, 100,0.,1. );
m_Energy_DEta[ieta] = new TH1F( TString("m_Energy_DEta")+index, TString("Two-cluster energy") +bin, 100,0.,20. );
m_Z_DEta[ieta] = new TH1F( TString("m_Z_DEta") +index, TString("Two cluster Z = |E#_1 - E#_2| / (E#_1+E#_2)")+bin, 100, 0.,1. );
mylist -> Add(m_Mass_DEta[ieta]);
mylist -> Add(m_Energy_DEta[ieta]);
mylist -> Add(m_Z_DEta[ieta]);
}
for ( Int_t iphi = 0; iphi < 60; iphi++ ) {
TString index = "["; index += iphi; index += "]";
TString bin = " #phi bin number "; bin += (iphi+1);
m_Mass_DPhi[iphi] = new TH1F( TString("m_Mass_DPhi") +index, TString("Two-cluster invariant mass")+bin, 100,0.,1. );
m_Energy_DPhi[iphi] = new TH1F( TString("m_Energy_DPhi")+index, TString("Two-cluster energy") +bin, 100,0.,20. );
m_Z_DPhi[iphi] = new TH1F( TString("m_Z_DPhi") +index, TString("Two cluster Z = |E#_1 - E#_2| / (E#_1+E#_2)")+bin, 100, 0.,1. );
mylist -> Add(m_Mass_DPhi[iphi]);
mylist -> Add(m_Energy_DPhi[iphi]);
mylist -> Add(m_Z_DPhi[iphi]);
}
for ( Int_t itow = 0; itow < 720; itow++ ) {
TString index = "["; index += itow; index += "]";
TString bin = " #tow bin number "; bin += (itow+1);
m_Mass_DTow[itow] = new TH1F( TString("m_Mass_DTow") +index, TString("Two-cluster invariant mass")+bin, 100,0.,1. );
m_Energy_DTow[itow] = new TH1F( TString("m_Energy_DTow")+index, TString("Two-cluster energy") +bin, 100,0.,20. );
m_Z_DTow[itow] = new TH1F( TString("m_Z_DTow") +index, TString("Two cluster Z = |E#_1 - E#_2| / (E#_1+E#_2)")+bin, 100, 0.,1. );
mylist -> Add(m_Mass_DTow[itow]);
mylist -> Add(m_Energy_DTow[itow]);
mylist -> Add(m_Z_DTow[itow]);
}
m_Preshower1_Adc1 = new TH1F("m_Preshower1_Adc1","Preshower-1 ADC spectrum for more energetic seed",4096,0.,4096.); // Preshower ADC spectrum, seed 1
m_Preshower2_Adc1 = new TH1F("m_Preshower2_Adc1","Preshower-2 ADC spectrum for more energetic seed",4096,0.,4096.); // Preshower ADC spectrum, seed 1
m_Postshower_Adc1 = new TH1F("m_Postshower_Adc1","Postshower ADC spectrum for more energetic seed",4096,0.,4096.); // Postshower ADC spectrum, seed 1
m_Preshower1_Adc1 = new TH1F("m_Preshower1_Adc2","Preshower-1 ADC spectrum for more energetic seed",4096,0.,4096.); // Preshower ADC spectrum, seed 2
m_Preshower2_Adc1 = new TH1F("m_Preshower2_Adc2","Preshower-2 ADC spectrum for more energetic seed",4096,0.,4096.); // Preshower ADC spectrum, seed 2
m_Postshower_Adc1 = new TH1F("m_Postshower_Adc2","Postshower ADC spectrum for more energetic seed",4096,0.,4096.); // Postshower ADC spectrum, seed 2
m_EndcapHits1 = new TH2F("m_EndcapHits1", "Location of higher-energy hit", 2*175, -175., 175., 2*175, -175., 175. );
m_EndcapHits2 = new TH2F("m_EndcapHits2", "Location of lower-energy hit", 2*175, -175., 175., 2*175, -175., 175. );
m_Energy1_vs_Energy2 = new TH2F("m_Energy1_vs_Energy2", "Energy of higher-energy hit vs lower-energy hit", 100, 0., 10., 100, 0., 10. );
m_NumClusters = new TH1F("m_NumClusters", "Number of clusters in the event", 100, 0., 100. );
m_NumSmdSeed15 = new TH1F("m_NumSmdSeed15","Number of SMD seeds, threshold 15 mips",6,0.,6.);
m_NumSmdSeed25 = new TH1F("m_NumSmdSeed25","Number of SMD seeds, threshold 25 mips",6,0.,6.);
m_NumSmdSeed35 = new TH1F("m_NumSmdSeed35","Number of SMD seeds, threshold 35 mips",6,0.,6.);
TIter next( mylist );
TH1 *h;
while ( (h = (TH1 *)next() ) ) {
h-> Sumw2();
}
}
/////////////////////////////////////////////////////////////////////////////
Int_t EEezPi0Analysis::Make( EEezClAnalysis *clusters,
EEstarTrig *trigger )
{
// Process an event using the cluster finder.
// Check trigger cuts
if ( !acceptEvent ( trigger ) ) return 1;
// Get a copy of out clusters
EEezClusterVec_t myClusters = clusters -> getClusters();
// Fill histogram w/ # of clusters
m_NumClusters -> Fill ( clusters -> getNClusters() );
// We're now ready to fill histograms...
Make ( &myClusters );
return 1;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
Int_t EEezPi0Analysis::Make ( EEezClusterVec_t *clusters ) {
//
// Iterate over all cluster-pairs, and calculate kinematic
// quantities assuming that the pair of clusters arises
// from pi0 --> gamma gamma.
//
// NOTE: This maker will assume that the cluster in the
// first position will always be cluster associated with
// the high-tower trigger. Ideally, this will be guarenteed
// by the cluster-maker.
// The first cluster in the list
EEezClusterVecIter_t iter1 = clusters -> begin();
// The maximum number of elements in the list to iterate
// over for the first iterator will be determined by the
// m_SearchType switch. We either search over all pairs
// of clusters, or we search over all pairs which include
// the tower associated with the high-tower cluster. In
// either case, the rest of the algorithmn remains the same.
EEezClusterVecIter_t lastCluster;
if ( m_SearchType == kAllPairs ) {
// First iterator will stop one short of the end
lastCluster = clusters -> end() - 1;
}
else if ( m_SearchType == kHighTowerOnly ) {
lastCluster = iter1 + 1;
}
else
assert(0); // Unknown type
// -----------<<<<<<<<<<< TOWER ONLY CLUSTER FINDER >>>>>>>>>>>-----------
// Nothing to do with towers only here.
if ( clusters -> size() < 2 ) return 1;
// Loop over all requested, unique cluster pairs in the list
do {
// Get the cluster (we do this for clarity...
// the compiler should optimize...)
EEezCluster cluster1 = (*iter1);
Float_t energy1 = cluster1 . getEnergy();
// The second iterator starts pointing one position
// further into the list than the current position
// of the first iterator, and continues to the end
// of the list.
EEezClusterVecIter_t iter2 = iter1 + 1;
do {
// The second cluster in the pair. NOTE:
// by construction, the first cluster _should_ have
// a more energetic seed tower than this cluster.
// this does not necessarily mean it has > energy
// (NEED TO HAVE CLUSTER FINDER ORDER CLUSTERS BY
// ADC, E AND/OR ET).
EEezCluster cluster2 = (*iter2);
Float_t energy2 = cluster2 . getEnergy();
if ( energy1 > energy2 )
FillHistograms( &cluster1, &cluster2 );
else
FillHistograms( &cluster2, &cluster1 );
} while ( ++iter2 != clusters -> end() );
} while ( ++iter1 != lastCluster );
// NOTE: The prefix-increment should be correct above, since
// we have just processed the nth element, and we want to make
// sure that the n+1th element is not the end.
return 1;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezPi0Analysis::FillHistograms ( EEezCluster *cluster1,
EEezCluster *cluster2,
Float_t weight )
{
// Fill histograms. Fill pi0 candidates.
//-- New Pi0 Candidate
EEezPi0nCandidate *candidate = new EEezPi0nCandidate();
candidate -> addCluster( cluster1 );
candidate -> addCluster( cluster2 );
Int_t nSmdSeed25 = candidate -> numberSmdSeeds(0,25.0);
nSmdSeed25 += candidate -> numberSmdSeeds(1,25.0);
m_NumSmdSeed25 -> Fill ( nSmdSeed25 );
Int_t nSmdSeed15 = candidate -> numberSmdSeeds(0,15.0);
nSmdSeed15 += candidate -> numberSmdSeeds(1,15.0);
m_NumSmdSeed15 -> Fill ( nSmdSeed15 );
Int_t nSmdSeed35 = candidate -> numberSmdSeeds(0,35.0);
nSmdSeed35 += candidate -> numberSmdSeeds(1,35.0);
m_NumSmdSeed35 -> Fill ( nSmdSeed35 );
//---
if ( nSmdSeed25 < mNSmdSeed ) return; // -- Require minimum of 3 SMD seeds
//
// Get the energy and momenta of the clusters
//
Float_t energy1 = cluster1 -> getEnergy();
TVector3 momentum1 = cluster1 -> getMomentum();
Int_t etabin1 = cluster1 -> getEtabin();
Int_t phibin1 = cluster1 -> getPhibin();
Int_t index1 = cluster1 -> getIndex();
Float_t energy2 = cluster2 -> getEnergy();
TVector3 momentum2 = cluster2 -> getMomentum();
Int_t etabin2 = cluster2 -> getEtabin();
Int_t phibin2 = cluster2 -> getPhibin();
Int_t index2 = cluster2 -> getIndex();
// Verify that the towers are within the allowed range
// for the calculation of kinematic quantities...
Int_t dEta = TMath::Abs( etabin1 - etabin2 );
Int_t dPhi = TMath::Min( TMath::Abs( phibin1 - phibin2 ), TMath::Abs( phibin1 - phibin2 - 60 ) );
// Veto pairs outside of patch m_Range on either side fo the
// more energetic tower...
if ( dEta > m_Range || dPhi > m_Range ) return;
// Veto pairs where seed towers are adjacent. Such
// combos will not be identified as separate clusters
// by the cluster finder.
if ( dEta <= 1 && dPhi <= 1 ) return;
// Calculate the relevant kinematic quantities
Float_t energy = energy1 + energy2;
assert( energy > 0. ); // or something is hosed
TVector3 momentum = momentum1 + momentum2;
Float_t mass = TMath::Sqrt( energy * energy - momentum * momentum );
Float_t Z = TMath::Abs( energy1 - energy2 ) / energy;
m_Mass -> Fill(mass, weight);
m_Mass1 -> Fill(mass, weight);
m_Mass2 -> Fill(mass, weight);
m_Mass_DEta[etabin1] -> Fill(mass, weight); m_Mass_DEta[etabin2] -> Fill(mass, weight);
m_Mass_DEta[phibin1] -> Fill(mass, weight); m_Mass_DEta[phibin2] -> Fill(mass, weight);
m_Mass_DEta[index1] -> Fill(mass, weight); m_Mass_DEta[index2] -> Fill(mass, weight);
//
// Fill quantities "triggered" on the proper mass range
//
if ( m_MinMass < mass && mass <= m_MaxMass ) {
m_Candidates.push_back(candidate);
m_Energy -> Fill(energy, weight);
m_Energy_DEta[etabin1] -> Fill(energy, weight); m_Energy_DEta[etabin2] -> Fill(energy, weight);
m_Energy_DEta[phibin1] -> Fill(energy, weight); m_Energy_DEta[phibin2] -> Fill(energy, weight);
m_Energy_DEta[index1] -> Fill(energy, weight); m_Energy_DEta[index2] -> Fill(energy, weight);
m_Z -> Fill(Z, weight);
m_Z_DEta[etabin1] -> Fill(Z, weight); m_Z_DEta[etabin2] -> Fill(Z, weight);
m_Z_DEta[phibin1] -> Fill(Z, weight); m_Z_DEta[phibin2] -> Fill(Z, weight);
m_Z_DEta[index1] -> Fill(Z, weight); m_Z_DEta[index2] -> Fill(Z, weight);
m_Eta -> Fill ( momentum.PseudoRapidity(), weight );
m_Phi -> Fill ( momentum.Phi(), weight );
m_Energy1_vs_Energy2 -> Fill ( energy2, energy1, weight );
Float_t x1 = ( kEEmcZSMD * cluster1 -> getMomentum().Unit() )[0];
Float_t y1 = ( kEEmcZSMD * cluster1 -> getMomentum().Unit() )[1];
Float_t x2 = ( kEEmcZSMD * cluster2 -> getMomentum().Unit() )[0];
Float_t y2 = ( kEEmcZSMD * cluster2 -> getMomentum().Unit() )[1];
m_EndcapHits1 -> Fill ( x1, y1 );
m_EndcapHits2 -> Fill ( x2, y2 );
}
else {
// Cleanup
delete candidate;
}
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
Int_t EEezPi0Analysis::Mix ( EEezClAnalysis *cluster_analy1,
EEstarTrig *trigger )
{
// Perform mixed-event background
// Check trigger cuts
if ( !acceptEvent( trigger ) ) return 1;
// Fill histogram w/ # of clusters
m_NumClusters -> Fill ( cluster_analy1 -> getNClusters() );
//
// Get the vectors of clusters from the cluster analysis
//
EEezClusterVec_t clusters1 = cluster_analy1 -> getClusters();
EEezClusterVec_t clusters2 = m_OldClusters;
if ( clusters1.size() == 0 ) return 1; //--
if ( clusters2.size() != 0 ) {
Float_t nc1 = (Float_t)clusters1.size();
Float_t nc2 = (Float_t)clusters2.size();
//
// Select a random cluster from the previous event
//
//$$$ Int_t irandom = (Int_t)(m_Random -> Rndm() * clusters2.size());
//
// Loop over all unique cluster pairs from the current event
//
EEezClusterVecIter_t iter1 = clusters1.begin();
while ( iter1 != clusters1.end()-1 ) {
EEezClusterVecIter_t iter2 = iter1+1;
while ( iter2 != clusters1.end() ) {
// Randomly select two clusters, one from this event
// and the other is the randomly selected cluster
// from the previous event.
Int_t iflip = ( m_Random -> Rndm() > 0.5 );
// Mix all clusters from previous event
for ( Int_t irandom = 0; irandom < (Int_t)clusters2.size(); irandom++ ) {
EEezCluster c1 = ( iflip ) ? (*iter1) : m_OldClusters[ irandom ];
EEezCluster c2 = ( !iflip ) ? (*iter2) : m_OldClusters[ irandom ];
if ( c1 . getEnergy() > c2 . getEnergy() )
FillHistograms( &c1, &c2, 1./nc2 );
else
FillHistograms( &c2, &c1, 1./nc2 );
} // all mix
iter2++;
}
iter1++;
}
}
// Copy clusters for next event
m_OldClusters.clear(); // we might get away without this
m_OldClusters = clusters1;
return 1;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
// Checks CTB sum and/or multiplicity
Int_t EEezPi0Analysis::acceptEvent( EEstarTrig *trigger )
{
// Checks CTB sum and/or multiplicity, and that the event
// satisfies the requested trigger.
// Return if the trigger-cut fails
if ( !trigger -> isTrigID( m_TriggerId ) ) return 0;
// Next verify that the CTB satisfies cuts.
Int_t n_ctb = 0;
Int_t sum_ctb = 0;
for ( Int_t i = 0; i < 240; i++ ) {
if ( (UShort_t)trigger -> CTB[i] > m_CtbThreshold ) {
n_ctb++;
m_CtbAdc -> Fill( (UShort_t)trigger->CTB[i] );
sum_ctb += trigger -> CTB[i];
}
}
// Cut on multiplicity
if ( n_ctb <= m_CtbMinMultiplicity ||
n_ctb > m_CtbMaxMultiplicity ) return 0;
// Cut on ADC sum
if ( sum_ctb <= m_CtbMinAdcSum ||
sum_ctb > m_CtbMaxAdcSum ) return 0;
// Return true
return 1;
};
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
void EEezPi0Analysis::Clear( Option_t *opts )
{
// Clear the vector of pi0n candidates
m_Candidates.clear();
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
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.