//////////////////////////////////////////////////////////////////////////
// //
// ATLFast MiscMaker class. //
// //
//////////////////////////////////////////////////////////////////////////
#include <TPythia.h>
#include <TMCParticle.h>
#include "ATLFMiscMaker.h"
#include "ATLFClusterMaker.h"
#include "ATLFPhotonMaker.h"
#include "ATLFMuonMaker.h"
#include "ATLFElectronMaker.h"
#include "ATLFJetMaker.h"
#include "ATLFMCMaker.h"
#include "ATLFCluster.h"
#include "ATLFElectron.h"
#include "ATLFJet.h"
#include "ATLFPhoton.h"
#include "ATLFMuon.h"
#include "ATLFMisc.h"
#include "ATLFast.h"
const Double_t kPI = TMath::Pi();
const Int_t nrludm = 100;
ClassImp(ATLFMiscMaker)
//_____________________________________________________________________________
ATLFMiscMaker::ATLFMiscMaker()
{
}
//_____________________________________________________________________________
ATLFMiscMaker::ATLFMiscMaker(const char *name, const char *title)
:ATLFMaker(name,title)
{
SetMinET();
m_Fruits = new ATLFMisc();
m_BranchName = "Misc";
Save();
}
//_____________________________________________________________________________
ATLFMiscMaker::~ATLFMiscMaker()
{
//dummy
}
//_____________________________________________________________________________
void ATLFMiscMaker::Clear(Option_t *option)
{
// Reset Misc Maker
ATLFMaker::Clear(option);
}
//_____________________________________________________________________________
void ATLFMiscMaker::Init()
{
// Create histograms
m_CircJets = new TH1F("CircJets","circularity of jets",50,0,1);
m_CircEvent = new TH1F("CircEvent","circularity of event",50,0,1);
m_Thrust = new TH1F("Thrust","thrust of event",50,0,1);
m_Oblateness = new TH1F("Oblateness","oblateness of event",50,0,1);
Int_t nbinc = 50;
Float_t tminc = -100;
Float_t tmaxc = 100;
m_RecPT = new TH1F("MisRecPT","reconstructed p_T",50,0,200);
m_RecPTcells = new TH1F("MisRecPTcells","reconstructed p_T +cells",50,0,200);
m_PTnu = new TH1F("MisPTnu","p_T nu",50,0,200);
m_PXnu = new TH1F("MisPXnu","pX nu",nbinc,tminc,tmaxc);
m_PYnu = new TH1F("MisPYnu","pY nu",nbinc,tminc,tmaxc);
m_RecPX = new TH1F("MisRecPX","reconstructed pX",nbinc,tminc,tmaxc);
m_RecPY = new TH1F("MisRecPY","reconstructed pY",nbinc,tminc,tmaxc);
m_RecPXcells = new TH1F("MisRecPXcells","reconstructed pX + cells",nbinc,tminc,tmaxc);
m_RecPYcells = new TH1F("MisRecPYcells","reconstructed pY + cells",nbinc,tminc,tmaxc);
m_ResolPX = new TH1F("MisResolPX","resol missing pX",nbinc,tminc,tmaxc);
m_ResolPX0 = new TH1F("MisResolPX0","resol missing pX (sumET=0-25)",nbinc,tminc,tmaxc);
m_ResolPX25 = new TH1F("MisResolPX25","resol missing pX (sumET=25-50)",nbinc,tminc,tmaxc);
m_ResolPX50 = new TH1F("MisResolPX50","resol missing pX (sumET=50-100)",nbinc,tminc,tmaxc);
m_ResolPX100 = new TH1F("MisResolPX100","resol missing pX (sumET=100-200)",nbinc,tminc,tmaxc);
m_ResolPX200 = new TH1F("MisResolPX200","resol missing pX (sumET=200-300)",nbinc,tminc,tmaxc);
m_ResolPX300 = new TH1F("MisResolPX300","resol missing pX (sumET=300-400)",nbinc,tminc,tmaxc);
Int_t nbina = 50;
Float_t tmina = 0.5;
Float_t tmaxa = 1.5;
m_PTmissPTnu = new TH1F("MisPTmissPTnu","PTmiss/PTnu",nbina,tmina,tmaxa);
m_PXmissPXnu = new TH1F("MisPXmissPXnu","PXmiss/PXnu",nbina,tmina,tmaxa);
m_PYmissPYnu = new TH1F("MisPYmissPYnu","PYmiss/PYnu",nbina,tmina,tmaxa);
}
//_____________________________________________________________________________
void ATLFMiscMaker::Finish()
{
// Function called by ATLFast::Finish at the end of the job
}
//_____________________________________________________________________________
Int_t ATLFMiscMaker::Make()
{
// Compute Miscelleneous info
MakeMisc();
// Compute missing energy
MakeMissing();
return 0;
}
//_____________________________________________________________________________
void ATLFMiscMaker::MakeMisc()
{
//.............................................
// This function calculates miscalleneous information about event itself:
// m_CircJets; //circularity of jets
// m_CircEvent; //circularity of event
// m_Thrust; //thrust of event
// m_Oblateness; //oblateness of event
//
// where m_CircJets and m_CircEvent denotes circularity calculated
// respectively from jets and from cells.
// m_Thrust and m_Oblateness represents respectively thrust axis of the event
// and oblateness.
//.............................................
ATLFMisc *misc = (ATLFMisc*)m_Fruits;
// Set Run and event number (taken from gATLFast object
misc->SetRun(gATLFast->Run());
misc->SetEvent(gATLFast->Event());
TPythia *pythia = (TPythia *)gATLFast->MCMaker()->Generator();
misc->SetMCProcess(pythia->GetMSTI(1));
// Set event characteristic
ATLFElectronMaker *electronmaker = gATLFast->ElectronMaker();
misc->SetNelectrons(electronmaker->Nelectrons());
ATLFPhotonMaker *photonmaker = gATLFast->PhotonMaker();
misc->SetNphotons(photonmaker->Nphotons());
ATLFJetMaker *jetmaker = gATLFast->JetMaker();
misc->SetNalljets(jetmaker->Njets());
TClonesArray *jets = jetmaker->Fruits();
ATLFJet *jet;
Int_t njets = jetmaker->Njets();
Int_t bjet = 0;
Int_t cjet = 0;
Int_t taujet = 0;
Int_t i;
for (i=0;i<njets;i++) {
jet = (ATLFJet*)jets->UncheckedAt(i);
if(TMath::Abs(jet->KFcode()) == 4) cjet++;
if(TMath::Abs(jet->KFcode()) == 5) bjet++;
if(TMath::Abs(jet->KFcode()) == 15) taujet++;
}
misc->SetNbjets(bjet);
misc->SetNcjets(cjet);
misc->SetNtaujets(taujet);
ATLFMuonMaker *muonmaker = gATLFast->MuonMaker();
TClonesArray *muons = muonmaker->Fruits();
ATLFMuon *muon;
Int_t nmuons = muonmaker->Nmuons();
Int_t muo = 0;
Int_t muox = 0;
for (i=0;i<nmuons;i++) {
muon = (ATLFMuon*)muons->UncheckedAt(i);
if(TMath::Abs(muon->Isolated()) == 1) muo++;
if(TMath::Abs(muon->Isolated()) != 1) muox++;
}
misc->SetNmuons(muo);
misc->SetNmuonsx(muox);
//........................................................
//..... Calculates circularity of jets and of event
//........................................................
Float_t circjet, circeve;
Circularity(circjet, circeve);
misc->SetCircJets(circjet);
misc->SetCircEvent(circeve);
m_CircJets->Fill(circjet);
m_CircEvent->Fill(circeve);
//........................................................
//.....calculates thrust and oblateness
//........................................................
Float_t thrust, obl;
ThrustOblateness(thrust, obl);
misc->SetThrust(thrust);
misc->SetOblateness(obl);
m_Thrust->Fill(thrust);
m_Oblateness->Fill(obl);
}
//_____________________________________________________________________________
void ATLFMiscMaker::MakeMissing()
{
//.............................................
// This function calculates the total missing transverse energy in the
// event. The transverse energies of cells not used for cluster
// reconstruction are smeared using the function ATLFClusterMaker::HadronResolution.
// The total measured transverse momenta p_xobserv and p_yobserv are obtained by
// summing over:
// jets, b-jets, c-jets, isolated leptons and photons, non-isolated muons which
// were not added to clusters, unused clusters, unused cells. The total missing
// transverse momenta pxmiss, pymiss are calculated as :
// pxmiss = - p_xobserv
// pymiss = - p_yobserv
// E_Tmiss = sqrt(pxmiss^2 + pymiss^2).
//
// The sum of the momenta of all particles escaping detection is
// also calculated. Results are stored in the object ATLFMisc, where
// in addition to pxmiss and pymiss above, pxnu and pynu
// denote the sum of momenta components of particles escaping detection
// (muons outside detector acceptance, neutrinos and SUSY LSP).
//.............................................
Int_t i;
//........................................................
//... Get pointers to MCparticles arrays and ClonesArray
//........................................................
ATLFMCMaker *mcarlo = gATLFast->MCMaker();
TClonesArray *particles = mcarlo->Fruits();
TMCParticle *part;
Int_t nparticles = particles->GetEntriesFast();
Float_t *pPT = mcarlo->PT();
Float_t *pEta = mcarlo->Eta();
//........................................................
//.... sum up reconstructed momenta
//........................................................
Float_t pxrec, pyrec, pxsum, pysum, sumet;
pxrec = pyrec = pxsum = pysum = sumet = 0;
//........................................................
//... Get pointers to jets arrays and ClonesArray
//...... add jets
//........................................................
ATLFJetMaker *jetmaker = gATLFast->JetMaker();
TClonesArray *jets = jetmaker->Fruits();
ATLFJet *jet;
Int_t njets = jetmaker->Njets();
for (i=0;i<njets;i++) {
jet = (ATLFJet*)jets->UncheckedAt(i);
pxrec += jet->PT()*TMath::Cos(jet->Phi());
pyrec += jet->PT()*TMath::Sin(jet->Phi());
sumet += jet->PT();
}
//........................................................
//... Get pointers to clustermaker arrays and ClonesArray
//...... add non-used clusters
//........................................................
ATLFClusterMaker *clustermaker = gATLFast->ClusterMaker();
TClonesArray *clusters = clustermaker->Fruits();
ATLFCluster *cluster;
Int_t nclusters = clustermaker->Nclusters();
Int_t ncells = clustermaker->Ncells();
for (i=0;i<nclusters;i++) {
cluster = (ATLFCluster*)clusters->UncheckedAt(i);
if (cluster->UseFlag() != 0) {
pxrec += cluster->ET()*TMath::Cos(cluster->Phi());
pyrec += cluster->ET()*TMath::Sin(cluster->Phi());
sumet += cluster->ET();
}
}
//........................................................
//... Get pointers to muonmaker arrays and ClonesArray
//.......add isolated muons
//........................................................
ATLFMuonMaker *muonmaker = gATLFast->MuonMaker();
TClonesArray *muons = muonmaker->Fruits();
ATLFMuon *muon;
Int_t nmuons = muonmaker->Nmuons();
for (i=0;i<nmuons;i++) {
muon = (ATLFMuon*)muons->UncheckedAt(i);
if (muon->Isolated() != 1) continue;
pxrec += muon->PT()*TMath::Cos(muon->Phi());
pyrec += muon->PT()*TMath::Sin(muon->Phi());
}
//........................................................
//.......add non-isolated muons not added to clusters
//........................................................
for (i=0;i<nmuons;i++) {
muon = (ATLFMuon*)muons->UncheckedAt(i);
if (muon->Isolated() != 0) continue;
if (muon->UseFlag() != 0 ) {
pxrec += muon->PT()*TMath::Cos(muon->Phi());
pyrec += muon->PT()*TMath::Sin(muon->Phi());
} else {
sumet += muon->PT();
}
}
//........................................................
//... Get pointers to electronmaker arrays and ClonesArray
//.......add isolated electrons
//........................................................
ATLFElectronMaker *electronmaker = gATLFast->ElectronMaker();
TClonesArray *electrons = electronmaker->Fruits();
ATLFElectron *electron;
Int_t nelectrons = electronmaker->Nelectrons();
for (i=0;i<nelectrons;i++) {
electron = (ATLFElectron*)electrons->UncheckedAt(i);
pxrec += electron->PT()*TMath::Cos(electron->Phi());
pyrec += electron->PT()*TMath::Sin(electron->Phi());
sumet += electron->PT();
}
//........................................................
//... Get pointers to photonmaker arrays and ClonesArray
//.......add isolated photons
//........................................................
ATLFPhotonMaker *photonmaker = gATLFast->PhotonMaker();
TClonesArray *photons = photonmaker->Fruits();
ATLFPhoton *photon;
Int_t nphotons = photonmaker->Nphotons();
for (i=0;i<nphotons;i++) {
photon = (ATLFPhoton*)photons->UncheckedAt(i);
pxrec += photon->PT()*TMath::Cos(photon->Phi());
pyrec += photon->PT()*TMath::Sin(photon->Phi());
sumet += photon->PT();
}
Float_t etrec = TMath::Sqrt(pxrec*pxrec + pyrec*pyrec);
m_RecPT->Fill(etrec);
m_RecPX->Fill(pxrec);
m_RecPY->Fill(pyrec);
//........................................................
//.......smear cells energy not used for reconstruction
//.......remove cells below threshold (after smearing)
//........................................................
Int_t *cellFlag = clustermaker->CellFlag();
Float_t *cellEta = clustermaker->CellEta();
Float_t *cellPhi = clustermaker->CellPhi();
Float_t *cellSumPT = clustermaker->CellSumPT();
Float_t pei, sigsme;
if (gATLFast->Smearing() == 1) {
for (i=0;i<ncells;i++) {
if (cellFlag[i] == 0) continue;
pei = cellSumPT[i]*TMath::CosH(cellEta[i]);
sigsme = clustermaker->HadronResolution(gATLFast->Luminosity(), pei, cellEta[i], cellSumPT[i],0);
pei *= 1+sigsme;
cellSumPT[i] = pei/TMath::CosH(cellEta[i]);
if (cellSumPT[i] < m_MinET) cellSumPT[i] = 0;
}
}
//........................................................
//.......add momenta of cells not used for reconstruction
//........................................................
for (i=0;i<ncells;i++) {
if (cellFlag[i] == 0) continue;
pxsum += cellSumPT[i]*TMath::Cos(cellPhi[i]);
pysum += cellSumPT[i]*TMath::Sin(cellPhi[i]);
sumet += cellSumPT[i];
}
pxsum += pxrec;
pysum += pyrec;
Float_t etsum = TMath::Sqrt(pxsum*pxsum + pysum*pysum);
m_RecPTcells->Fill(etsum);
m_RecPXcells->Fill(pxsum);
m_RecPYcells->Fill(pysum);
Float_t pxmiss = -pxsum;
Float_t pymiss = -pysum;
Float_t ptmiss = TMath::Sqrt(pxmiss*pxmiss + pymiss*pymiss);
//........................................................
//......sum up momenta of neutrinos + SUSY LSP
//........................................................
Float_t pxnu = 0;
Float_t pynu = 0;
Int_t akf;
Int_t kFLSP = gATLFast->SUSYcodeLSP();
for (i=mcarlo->Nstart();i<nparticles;i++) {
part = (TMCParticle*)particles->UncheckedAt(i);
akf = TMath::Abs(part->GetKF());
if (akf == 12 || akf == 14 || akf == 16 || akf == kFLSP) {
pxnu += part->GetPx();
pynu += part->GetPy();
}
//........................................................
//......add muons which are outside acceptance to the missing momenta
//........................................................
if (akf == 13) {
if (pPT[i] < muonmaker->MinPT() || TMath::Abs(pEta[i]) > muonmaker->MaxEta()) {
pxnu += part->GetPx();
pynu += part->GetPy();
}
}
}
Float_t ptnu = TMath::Sqrt(pxnu*pxnu + pynu*pynu);
m_PTnu->Fill(ptnu);
m_PXnu->Fill(pxnu);
m_PYnu->Fill(pynu);
if (ptnu > 0) m_PTmissPTnu->Fill(ptmiss/ptnu);
if (pxnu > 0) m_PXmissPXnu->Fill(ptmiss/pxnu);
if (pynu > 0) m_PYmissPYnu->Fill(ptmiss/pynu);
//........................................................
//.......monitoring p_T_miss resolution in sumE_T bins
//........................................................
if (sumet > 0 && sumet < 25) m_ResolPX0->Fill(pxmiss-pxnu);
if (sumet > 25 && sumet < 50) m_ResolPX25->Fill(pxmiss-pxnu);
if (sumet > 50 && sumet < 100) m_ResolPX50->Fill(pxmiss-pxnu);
if (sumet > 100 && sumet < 200) m_ResolPX100->Fill(pxmiss-pxnu);
if (sumet > 200 && sumet < 300) m_ResolPX200->Fill(pxmiss-pxnu);
if (sumet > 300 && sumet < 400) m_ResolPX300->Fill(pxmiss-pxnu);
m_ResolPX->Fill(pxmiss-pxnu);
//........................................................
//.......Store missing energy parameters in Root object
//........................................................
ATLFMisc *misc = (ATLFMisc*)m_Fruits;
misc->SetMissing(pxmiss,pymiss,pxnu,pynu);
}
//_____________________________________________________________________________
void ATLFMiscMaker::PrintInfo()
{
ATLFMaker::PrintInfo();
}
//_____________________________________________________________________________
void ATLFMiscMaker::Circularity(Float_t &circj, Float_t &circev)
{
//........................................................
//...Purpose: to calculate sphericity of the event
// Original routine SPHERIC from LUJET (and E. Richter-Was Spring 1996)
//
// J Soderqvist 09/12/96
// Rewrite to calculate Circularity of all jets (CIRCJ)
// and of event (CIRCEV).
// Use gen routine eisrs1.f to calculate eigenvalues
//
//........................................................
//------ matrices to be filled
Float_t circ1[2][2], circ2[2][2];
//------ eigenvalues and eigenvectors to be calculated
Float_t circval1[2], circval2[2], circvec1[2][2], circvec2[2][2];
Float_t work[2];
Int_t i, j, j1,j2;
Float_t p[2];
//------ zero all vectors and matrices
for (i=0;i<2;i++) {
circval1[i] = 0;
circval2[i] = 0;
for (j=0;j<2;j++) {
circ1[i][j] = 0;
circ2[i][j] = 0;
circvec1[i][j] = 0;
circvec2[i][j] = 0;
}
}
//------ Calculate matrix to be diagonalized for all calorimeter
ATLFClusterMaker *clustermaker = gATLFast->ClusterMaker();
Int_t ncells = clustermaker->Ncells();
Float_t *cellPhi = clustermaker->CellPhi();
Float_t *cellSumPT = clustermaker->CellSumPT();
for (i=0;i<ncells;i++) {
p[0] = cellSumPT[i]*TMath::Cos(cellPhi[i]);
p[1] = cellSumPT[i]*TMath::Sin(cellPhi[i]);
for (j1=0;j1<2;j1++) {
for (j2=0;j2<2;j2++) {
circ1[j1][j2] += p[j1]*p[j2];
}
}
}
//------ Calculate Circularity for CELLS in Calorimeter
Int_t ierr = 0;
EigenValues(2,2,circ1,circval1,circvec1,ierr,work);
if (ierr) {
Warning("Circularity","in call to eisrs1, CIRC1");
circev = -777;
return;
}
if (circval1[0]+circval1[1] != 0) {
circev = 2*TMath::Min(circval1[0], circval1[1])/(circval1[0] + circval1[1]);
}
//------ Calculate matrix to be diagonalized for JETS
ATLFJetMaker *jetmaker = gATLFast->JetMaker();
TClonesArray *jets = jetmaker->Fruits();
ATLFJet *jet;
Int_t njets = jetmaker->Njets();
//....include jets
for (i=0;i<njets;i++) {
jet = (ATLFJet*)jets->UncheckedAt(i);
p[0] = jet->PT()*TMath::Cos(jet->Phi());
p[1] = jet->PT()*TMath::Sin(jet->Phi());
for (j1=0;j1<2;j1++) {
for (j2=0;j2<2;j2++) {
circ2[j1][j2] += p[j1]*p[j2];
}
}
}
//------ Calculate Circularity for JETS
EigenValues(2,2,circ2,circval2,circvec2,ierr,work);
if (ierr) {
Warning("Circularity","in call to eisrs1, CIRC2");
circj = -777;
return;
}
if (circval2[0]+circval2[1] != 0) {
circj = 2*TMath::Min(circval2[0], circval2[1])/(circval2[0] + circval2[1]);
}
//printf("Exiting Circularity, cirev=%g, circj=%gn",circev,circj);
}
//_____________________________________________________________________________
void ATLFMiscMaker::ThrustOblateness(Float_t &thrust, Float_t &obl)
{
//........................................................
// J Soderqvist 10/12/96
//
// Function to calculate thrust and oblateness of events
// in ATLFAST. Take clusters from ATLFAST and use
// code that M. Pearce originally extracted from OPAL
// to calculate thrust and oblateness.
// (ATLUT3 is using JETSET thrust finding algorithm,
// but is not dependent on the JETSET commons)
//
// Modified 30/12/96 J.S.
// to give values -777 if NCLU=0
//
// The parameters to pass to ATLUT3 are:
//
// function ATLUT3 (N,NRLUDM,P,THR,OBL)
//
// INPUT : n Number of tracks
// INPUT : nrludm First dimension of PLUND
// INPUT : p 4-momenta in Jetset format (Px,Py,Pz,E)
// OUTPUT : thr Thrust value
// OUTPUT : obl Oblateness value
//
// Use clusters identified by ATLFAST
//........................................................
Int_t i,k;
Float_t patl[nrludm][5], thratl, oblatl;
ATLFClusterMaker *clustermaker = gATLFast->ClusterMaker();
TClonesArray *clusters = clustermaker->Fruits();
ATLFCluster *cluster;
Int_t nclusters = clustermaker->Nclusters();
for (i=0;i<nrludm;i++) {
for (k=0;k<5;k++) patl[i][k] = 0;
}
if (nclusters > 0) {
for (i=0;i<nclusters;i++) {
cluster = (ATLFCluster*)clusters->UncheckedAt(i);
patl[i][0] = cluster->ET()*TMath::Cos(cluster->Phi());
patl[i][1] = cluster->ET()*TMath::Sin(cluster->Phi());
patl[i][2] = cluster->ET()*TMath::SinH(cluster->Eta());
patl[i][3] = TMath::Sqrt(patl[0][i]*patl[i][0] + patl[i][1]*patl[i][1] + patl[i][2]*patl[i][2]);
patl[i][4] = 0;
}
thratl = -1;
oblatl = -1;
// Atlut3(ncluatl, patl, thratl, oblatl);
//printf("ThrustOblateness: nclusters=%d, thratl=%g, oblatl=%gn",nclusters,thratl,oblatl);
thrust = thratl;
obl = oblatl;
} else {
thrust = -777;
obl = -777;
}
}
//#ifdef NEWCODE
//_____________________________________________________________________________
void ATLFMiscMaker::EigenValues(Int_t nm, Int_t n, Float_t ar[2][2], Float_t *wr, Float_t zr[2][2],
Int_t &ierr, Float_t *work)
{
// Compute all eigenvalues and corresponding eigenvectors of a real symmetric matrix.
// translated to C++ from CERNLIB routine EISRS1
EigenValue1(nm,n,ar,wr,work,zr);
EigenValue2(nm,n,wr,work,zr,ierr);
}
//_____________________________________________________________________________
void ATLFMiscMaker::EigenValue1(Int_t , Int_t n, Float_t a[2][2], Float_t *d, Float_t *e, Float_t z[2][2])
{
// translated to C++ from CERNLIB routine TRED2
Int_t i,j,k,l,ii,jp1;
Float_t f,g,h,hh,scale;
for (i=0;i<n;i++) {
for (j=0;j<i;j++) z[i][j] = a[i][j];
}
if (n == 1) goto L320;
for (ii=1;ii<n;ii++) {
i = n+1-ii;
l = i-1;
h = 0;
scale = 0;
if (l < 1) {e[i] = z[i][l]; d[i] = h; continue;}
for (k=0;k<l;k++) scale += TMath::Abs(z[i][k]);
if (scale == 0) { e[i] = z[i][l]; d[i] = h; continue;}
for (k=0;k<l;k++) {
z[i][k] /= scale;
h += z[i][k]*z[i][k];
}
f = z[i][l];
g = -TMath::Sign(Float_t(TMath::Sqrt(h)), f);
e[i] = scale*g;
h -= f*g;
z[i][l] = f-g;
f = 0;
for (j=0;j<l;j++) {
z[j][i] /= scale*h;
g = 0;
for (k=0;k<j;k++) g += z[j][k]*z[i][k];
jp1 = j+1;
if (l >= jp1) {
for (k=jp1;k<l;k++) g += z[k][j]*z[i][k];
}
e[j] = g/h;
f += e[j]*z[i][j];
}
hh = f/(h+h);
for (j=0;j<l;j++) {
f = z[i][j];
g = e[j] -hh*f;
e[j] = g;
for (k=0;k<j;k++) z[i][k] += -f*e[k] -g*z[i][k];
}
for (k=0;k<l;k++) z[i][k] *= scale;
d[i] = h;
}
L320:
d[0] = 0;
e[0] = 0;
for (i=0;i<n;i++) {
l = i-1;
if (d[i] != 0) {
for (j=0;j<l;j++) {
g = 0;
for (k=0;k<l;k++) g += z[i][k]*z[k][j];
for (k=0;k<l;k++) z[k][j] += -g*z[k][i];
}
}
d[i] = z[i][i];
z[i][i] = 1;
if (l < 0) continue;
for (j=0;j<l;j++) { z[i][j] = 0; z[j][i] = 0; }
}
}
//_____________________________________________________________________________
void ATLFMiscMaker::EigenValue2(Int_t , Int_t n, Float_t *d, Float_t *e, Float_t z[2][2], Int_t &ierr)
{
// translated to C++ from CERNLIB routine TQL2
Int_t i,j,k,l,m,ii,mml;
Float_t b,c,f,g,h,p,r,s;
const Float_t MACHEP = TMath::Power(2,-23);
ierr = 0;
if (n == 1) return;
for (i=1;i<n;i++) e[i-1] = e[i];
f = b = 0;
e[n-1] = 0;
for (l=0;l<n;l++) {
j = 0;
h = MACHEP*(TMath::Abs(d[l]) + TMath::Abs(e[l]));
if (b < h) b = h;
for (m=l;m<n;m++) {
if (TMath::Abs(e[m]) <= b) break;
}
if (m == l) {d[l] += f; continue;}
while(1) {
if (j == 30) {ierr = l; return;}
j++;
p = (d[l+1]-d[l])/(2*e[l]);
r = TMath::Sqrt(p+p+1);
h = d[l] -e[l]/(p+TMath::Sign(r,p));
for (i=l;i<n;i++) d[i] -= h;
f += h;
p = d[m];
c = 1;
s = 0;
mml = m-l;
for (ii=0;ii<mml;ii++) {
i = m-ii;
g = c*e[i];
h = c*p;
if (TMath::Abs(p) < TMath::Abs(e[i])) {
c = p/e[i];
r = TMath::Sqrt(c*c+1);
e[i+1] = s*e[i]*r;
s = 1/r;
c *= s;
} else {
c = e[i]/p;
r = TMath::Sqrt(c*c+1);
e[i+1] = s*p*r;
s = c/r;
c = 1/r;
}
p = c*d[i] - s*g;
d[i+1] = h + s*(c*g + s*d[i]);
for (k=0;k<n;k++) {
h = z[k][i+1];
z[k][i+1] = s*z[k][i] +c*h;
z[k][i] = c*z[k][i] - s*h;
}
}
e[l] = s*p;
d[l] = c*p;
if (TMath::Abs(e[l]) <= b) {d[l] += f; break;}
}
}
for (ii=1;ii<n;ii++) {
i = ii-1;
k = i;
p = d[i];
for (j=ii;j<n;j++) {
if (d[j] >= p) continue;
k = j;
p = d[j];
}
if (k == i) continue;
d[k] = d[i];
d[i] = p;
for (j=0;j<n;j++) {
p = z[j][i];
z[j][i] = z[j][k];
z[j][k] = p;
}
}
}
//_____________________________________________________________________________
void ATLFMiscMaker::Atlut3(Int_t n,Float_t p[][5], Float_t &thr, Float_t &obl)
{
//........................................................
//. Extracted by M. Pearce from OPAL physics utilities (PX library)
//. Originally authored by W. J. Gary
//. This version converted for ATLAS use by
//. M. Pearce and J. Soderqvist, KTH Stockholm, December 1996
//.
//. ------
//. ATLUT3
//. ------
//. An "in-house" version of the Jetset thrust finding algorithm
//. which works entirely through an argument list rather than
//. with e.g. the Jetset common blocks. This routine calculates
//. the standard linear thrust vectors and values. Its operation
//. is entirely decoupled from any MST or MSTJ variables etc.
//. which might be set by a user using Jetset.
//. The main purpose of developing an in-house version of the
//. Jetset thrust algorithm is so as to have a version
//. which is compatible with both Jetset6.3 and Jetset7.1 etc.
//. (because of the change in the Jetset common blocks between
//. these two versions, the Jetset version of this thrust
//. algorithm LUTHRU is version specific).
//.
//. The Jetset thrust algorithm implements an "approximate" method
//. for thrust determination because not all particle combinations
//. are tested. It is therefore logically possible that the thrust
//. axes and values determined by this routine will correspond
//. to local rather than to absolute maxima of the thrust function.
//. However in practice this is unlikely because several starting
//. points are used and the algorithm iterated to cross check one
//. convergence vs. another. Thus this routine offers a very fast
//. and in practice quite accurate algorithm for thrust (much faster
//. than so-called "exact" algorithms).
//. Usage :
//.
//.
//. INPUT : ntrak Number of tracks
//. INPUT : nrludm First dimension of PLUND
//. INPUT : plund 4-momenta in Jetset format
//. OUTPUT : thrust Thrust value
//. OUTPUT : oblate Oblateness value
//.
//. CALLS : Atanxy, Atplu3, Atrmx3, Atrof3, Atrob3
//. CALLED : Atlth4
//.
//. AUTHOR : Modified from LUTHRU (T.Sjostrand) by J.W.Gary (OPAL)
//. CREATED : 31-Jan-89
//. LAST MOD : Dec-96
//.
//. OPAL Modification Log.
//. 04-Feb-89 In-house version for PX library J.W.Gary
//. 12-Mar-89 Get rid of calls to RLU J.W.Gary
//. 27-Nov-95 M.Schroder Clear part of the array P above tracks
//. ATLAS Modification Log
//. Dec-96 ATLAS naming of routines M. Pearce and J. Soderqvist
//........................................................
Int_t np,ilc,ild,ilf,ilg,i,j,nc,ipp,ierr;
Float_t tdi[3], tpr[3], pvec[3], rvec[3], rmx[3][3];
Float_t ps,tds,sgn,thp,thps,phi,cp,sp;
Int_t iagr = 0;
const Float_t paru42 = 1;
const Float_t paru48 = 0.0001;
const Int_t mstu44 = 4;
const Int_t mstu45 = 2;
if (2*n+mstu44+15 >= nrludm-5) {
thr = -2;
obl = -2;
return;
}
// M.Schroder (these elements are always used, but sometimes not set...)
// for (i=n;i<2*n+2;i++) {
// p[i][0] = 0;
// p[i][1] = 0;
// p[i][2] = 0; // <===already preset in calling routine
// p[i][3] = 0;
// p[i][4] = 0;
// }
//...Take copy of particles that are to be considered in thrust analysis.
np = 0;
ps = 0;
for (i=0;i<n;i++) {
p[n+np][0] = p[i][0];
p[n+np][1] = p[i][1];
p[n+np][2] = p[i][2];
p[n+np][3] = TMath::Sqrt(p[i][0]*p[i][0] +p[i][1]*p[i][1] + p[i][2]*p[i][2]);
p[n+np][4] = 1;
if (TMath::Abs(paru42-1) > 0.001) p[n+np][4] = TMath::Power(p[n+np][3], paru42-1);
ps += p[n+np][3]*p[n+np][4];
np++;
}
//...Loop over thrust and major. T axis along z direction in latter case.
for (ild=0;ild<2;ild++) {
if (ild == 1) {
Atanxy(p[n+np][0], p[n+np][1], phi, ierr);
Atplu3(n+np+1, p, pvec, "U");
cp = TMath::Cos(phi);
sp = TMath::Sin(phi);
Atrmx3(pvec, cp, sp, rmx);
for (ipp=n;ipp<n+np+1;ipp++) {
Atplu3(ipp,p,pvec,"U");
Atrof3(rmx,pvec,rvec);
Atplu3(ipp,p,rvec,"P");
}
}
//...Find and order particles with highest p (pT for major).
for (ilf=n+np+3;ilf<n+np+mstu44+4;ilf++) p[ilf][3] = 0;
for (i=n;i<n+np;i++) {
if (ild == 1) p[i][3] = TMath::Sqrt(p[i][0]*p[i][0] + p[i][1]*p[i][1]);
for (ilf=n+np+mstu44+2;ilf>=n+np+4;ilf--) {
if (p[i][3] <= p[ilf][3]) goto L130;
for (j=0;j<5;j++) p[ilf+1][j] = p[ilf][j];
}
ilf = n+np+3;
L130:
for(j=0;j<5;j++) p[ilf][j] = p[i][j];
}
//...Find and order initial axes with highest thrust (major).
for (ilg=n+np+mstu44+4;ilg<n+np+mstu44+15;ilg++) p[ilg][3] = 0;
nc = Int_t(TMath::Power(2, TMath::Min(mstu44,np)-1));
Int_t ilf1, ilf2;
for (ilc=0;ilc<nc;ilc++) {
for (j=0;j<3;j++) tdi[j] = 0;
for (ilf=0;ilf<TMath::Min(mstu44,np);ilf++) {
sgn = p[n+np+ilf+3][4];
ilf1 = Int_t(TMath::Power(2,ilf));
ilf2 = 2*ilf1;
if (ilf2*((ilc+ilf1-1)/ilf2) >= ilc+1) sgn = -sgn;
for (j=0;j<3-ild;j++) tdi[j] += sgn*p[n+np+ilf+3][j];
}
tds = tdi[0]*tdi[0] + tdi[1]*tdi[1] + tdi[2]*tdi[2];
for (ilg=n+np+mstu44+TMath::Min(ilc,9)+4;ilg<n+np+mstu44+5;ilg--) {
if (tds <= p[ilg][3]) goto L200;
for (j=0;j<4;j++) p[ilg+1][j] = p[ilg][j];
}
ilg = n+np+mstu44 +4;
L200:
for (j=0;j<3;j++) p[ilg][j] = tdi[j];
p[ilg][3] = tds;
}
//...Iterate direction of axis until stable maximum.
p[n+np+ild][3] = 0;
ilg = 0;
while(1) {
ilg++;
thp = 0;
while(1) {
thps = thp;
for (j=0;j<3;j++) {
if (thp <= 1e-10) tdi[j] = p[n+np+mstu44+3+ilg][j];
else tdi[j] = tpr[j];
tpr[j] = 0;
}
for (i=n;i<n+np;i++) {
sgn = TMath::Sign(p[i][4], tdi[0]*p[i][0] +tdi[1]*p[i][1] + tdi[2]*p[i][2]);
for (j=0;j<3-ild;j++) tpr[j] += sgn*p[i][j];
}
thp = TMath::Sqrt(tpr[0]*tpr[0] + tpr[1]*tpr[1] + tpr[2]*tpr[2]);
if (thp < thps+paru48) break;
}
//...Save good axis. Try new initial axis until a number of tries agree.
if (thp < p[n+np+ild][3]-paru48 && ilg < TMath::Min(10,nc)) continue;
if (thp > p[n+np+ild][3]+paru48) {
iagr = 0;
sgn = 1;
for (j=0;j<3;j++) p[n+np+ild][j] = sgn*tpr[j]/(ps*thp);
p[n+np+ild][3] = thp;
p[n+np+ild][4] = 0;
}
iagr++;
if (iagr >= mstu45 || ilg >= TMath::Min(10,nc)) break;
}
}
//...Find minor axis and value by orthogonality.
//**JWG SGN = (-1.)**INT (RLU(0)+0.5)
sgn = 1;
p[n+np+2][0] = -sgn*p[n+np+1][1];
p[n+np+2][1] = sgn*p[n+np+1][0];
p[n+np+2][2] = 0;
thp = 0;
for (i=n;i<n+np;i++) {
thp += p[i][4]*TMath::Abs(p[n+np+2][0]*p[i][0] + p[n+np+2][1]*p[i][1]);
}
p[n+np+2][3] = thp/ps;
p[n+np+2][4] = 0;
//...Fill axis information. Rotate back to original coordinate system.
for (ild=0;ild<3;ild++) {
for (j=0;j<5;j++) p[n+ild][j] = p[n+np+ild][j];
}
for (ipp=n;ipp<n+3;ipp++) {
Atplu3(ipp,p,pvec,"U");
Atrob3(rmx,pvec,rvec);
Atplu3(ipp,p,rvec,"P");
}
//...Select storing option. Calculate thurst and oblateness.
thr = p[n][3];
obl = p[n+1][3] - p[n+2][3];
}
//_____________________________________________________________________________
void ATLFMiscMaker::Atanxy(Float_t xx, Float_t yy, Float_t &ang, Int_t &ierr)
{
//........................................................
//. SOURCE: Jetset7.1 (T. Sjostrand)
//. Reconstruct the azimuthal angle of a vector,
//. given the X and Y components of a vector
//. Usage :
//.
//.
//. INPUT : xx The X component of a vector
//. INPUT : yy The Y component of a vector
//. OUTPUT : Ang The azimuthal angle
//. OUTPUT : ierr = 0 if all is OK ; = -1 otherwise
//........................................................
Double_t ulangl, rrr, xxx, yyy;
ierr = 0;
xxx = xx;
yyy = yy;
rrr = TMath::Sqrt(xxx*xxx + yyy*yyy);
if (rrr < 1e-20) {ierr = -1; return;}
if (TMath::Abs(xxx)/rrr < 0.8) {
ulangl = TMath::Sign(TMath::ACos(xxx/rrr), yyy);
} else {
ulangl = TMath::ASin(yyy/rrr);
if (xxx < 0 && ulangl >= 0) ulangl = kPI - ulangl;
else if (xxx < 0) ulangl = -kPI -ulangl;
}
ang = ulangl;
}
#ifdef NEW
//_____________________________________________________________________________
void ATLFMiscMaker::Atlth4(Int_t ntrak, Int_t itkdm, Float_t *ptrak, Float_t *thrval, Float_t thrvec[][3], Int_t &ierr)
{
//........................................................
//. Extracted by M. Pearce from OPAL physics utilities (PX library)
//. Originally authored by W. J. Gary
//. This version converted for ATLAS use by
//. M. Pearce and J. Soderqvist, KTH Stockholm, December 1996
//.
//. ------
//. ATLTH4
//. ------
//. Routine to determine the Thrust Principal, Major and
//. Minor axes and values using the Jetset algorithm.
//. The implementation here is without a common block, however.
//. Thus this routine may be used regardless of whether the
//. Jetset6.3 or Jetset7.1 library might be linked. It is
//. not necessary to link to Jetset, however.
//. Usage :
//.
//. INTEGER ITKDM,MXTRAK
//. PARAMETER (ITKDM=3.or.more,MXTRAK=1.or.more)
//. INTEGER NTRAK,IERR
//. REAL PTRAK (ITKDM,MXTRAK),
//. + THRVEC (3,3.or.more),
//. + THRVAL (3.or.more)
//.
//. NTRAK = 1.to.MXTRAK
//. CALL ATLTH4 (NTRAK,ITKDM,PTRAK,THRVAL,THRVEC,IERR)
//.
//. The thrust vectors THRVEC are ordered according to the
//. corresponding thrust values THRVAL such that
//. THRVAL (1) < THRVAL (2) < THRVAL (3)
//. Thus THRVEC (*,3) is the Thrust Principal axis;
//. Thus THRVEC (*,2) is the Thrust Major axis;
//. Thus THRVEC (*,1) is the Thrust Minor axis;
//.
//. INPUT : ntrak Total number of particles
//. INPUT : itkdm First dimension of PTRAK array
//. INPUT : ptrak Particle momentum array: Px,Py,Pz,E
//. OUTPUT : thrval Thrust values
//. OUTPUT : thrvec Associated Thrust vectors
//. OUTPUT : ierr = 0 if all is OK ; = -1 otherwise
//.
//. AUTHOR : J.W.Gary
//. CREATED : 18-Jun-88
//. LAST MOD : Dec-96
//.
//. OPAL Modification Log.
//. 04-Feb-89 Integrate with ATLUT3 J.W.Gary
//. ATLAS Modification Log.
//. Dec-96 ATLAS naming of routines M. Pearce and J. Soderqvist
//........................................................
Int_t axis,ix1, ix2;
Float_t thrust, oblate;
ierr = 0;
if (ntrak <= 1 || ntrack > nrludm) {ierr = -1; return;}
//* Pack 4-momenta in Jetset format
//* ---- --------- -- ------ ------
for (ix1=0;ix1<ntrak;ix1++) {
for (axis=0;axis<4;axis++) {
plund[ix1][axis] = ptrak[axis][ix1];
}
}
//* Jetset algorithm for Thrust
//* ------ --------- --- ------
Atlut3(ntrak,plund,thrust,oblate);
if (thrust < 0) {ierr = -1; return;}
//* Buffer eigenvectors for output
//* ------ ------------ --- ------
for (ix1=0;ix1<3;ix1++) {
ix2 = ntrak + (3-ix1) -1;
thrval[ix1] = plund[ix2][3];
for (axis=0;axis<4;axis++) {
thrvec[axis][ix1] = plund[ix2][axis];
}
}
}
#endif
//_____________________________________________________________________________
void ATLFMiscMaker::Atplu3(Int_t indx, Float_t plund[][5], Float_t *pvec, char *option)
{
//........................................................
//. Extracted by M. Pearce from OPAL physics utilities (PX library)
//. Originally authored by W. J. Gary
//. This version converted for ATLAS use by
//. M. Pearce and J. Soderqvist, KTH Stockholm, December 1996
//.
//. ------
//. ATPLU3
//. ------
//. A utility routine to repack a Jetset 3-momentum as a
//. standard ordered array or vice-versa. This routine
//. is used to translate between the array conventions
//. of the Jetset thrust algorithm and of the other routines
//. in this library
//. Usage :
//.
//.
//. INPUT : indx The Jetset vector number
//. IN/OUT : nrludm First argument of PLUND
//. IN/OUT : plund The Jetset 3-momentum array
//. IN/OUT : pvec The input or output array
//. CONTROL : option ='U' means unpack Jetset array
//. = anything else means pack Jetset array
//.
//. AUTHOR : J.W.Gary
//. CREATED : 04-Feb-89
//. LAST MOD : 04-Feb-89
//.
//. ATLAS Modification Log.
//. Dec-96 ATLAS naming of routines M. Pearce and J. Soderqvist
//........................................................
Int_t ix;
for (ix=0;ix<3;ix++) {
if (option[0] == 'U') pvec[ix] = plund[indx][ix];
else plund[indx][ix] = pvec[ix];
}
}
//_____________________________________________________________________________
void ATLFMiscMaker::Atrob3(Float_t rmx[3][3], Float_t *vect, Float_t *rvec)
{
//........................................................
//. SOURCE: HERWIG (B.Webber,G.Marchesini)
//. Rotate 3-vector VECT by inverse of rotation matrix RMX,
//. RVEC = (RMX)-1 * VECT
//. Usage :
//.
//. Float_t vect (3.or.more),
//. + rvec (3.or.more),
//. + rmx (3,3.or.more)
//.
//.
//. INPUT : rmx The rotation matrix
//. INPUT : vect The vector to be rotated
//. OUTPUT : rvec The rotated vector
//........................................................
Double_t s1, s2, s3;
s1 = vect[0]*rmx[0][0] + vect[1]*rmx[1][0] + vect[2]*rmx[2][0];
s2 = vect[0]*rmx[0][1] + vect[1]*rmx[1][1] + vect[2]*rmx[2][1];
s3 = vect[0]*rmx[0][2] + vect[1]*rmx[1][2] + vect[2]*rmx[2][2];
rvec[0] = s1;
rvec[1] = s2;
rvec[2] = s3;
}
//_____________________________________________________________________________
void ATLFMiscMaker::Atrof3(Float_t rmx[3][3], Float_t *vect, Float_t *rvec)
{
//........................................................
//. SOURCE: HERWIG (B.Webber,G.Marchesini)
//. Rotate 3-vector VECT by rotation matrix RMX,
//. RVEC = RMX * VECT
//. Usage :
//.
//. Float_t vect (3.or.more),
//. + rvec (3.or.more),
//. + rmx (3,3.or.more)
//.
//. CALL ATROF3 (RMX,VECT,RVEC)
//.
//. INPUT : rmx The rotation matrix
//. INPUT : vect The vector to be rotated
//. OUTPUT : rvec The rotated vector
//........................................................
Double_t s1, s2, s3;
s1 = rmx[0][0]*vect[0] + rmx[0][1]*vect[1] + rmx[0][2]*vect[2];
s2 = rmx[1][0]*vect[0] + rmx[1][1]*vect[1] + rmx[1][2]*vect[2];
s3 = rmx[2][0]*vect[0] + rmx[2][1]*vect[1] + rmx[2][2]*vect[2];
rvec[0] = s1;
rvec[1] = s2;
rvec[2] = s3;
}
//_____________________________________________________________________________
void ATLFMiscMaker::Atrmx3(Float_t *vect, Float_t cp, Float_t sp, Float_t rmx[3][3])
{
//........................................................
//. SOURCE: HERWIG (B.Webber,G.Marchesini)
//. Calculate the rotation matrix to get from vector VECT
//. to the Z axis, followed by a rotation PHI about the Z axis,
//. where CP, SP are the cosine and sine of PHI, respectively.
//. Usage :
//.
//. INPUT : vect The vector for which the rotation matrix
//. is to be calculated
//. INPUT : cp Cosine of phi
//. INPUT : sp Sine of phi
//. OUTPUT : rmx The rotation matrix
//........................................................
Double_t ct, st, cf, sf, pp, pt;
const Float_t ptcut = 1e-10;
pt = vect[0]*vect[0] + vect[1]*vect[1];
if (pt < ptcut) {
ct = TMath::Sign(Float_t(1), vect[2]);
st = 0;
cf = 1;
sf = 0;
} else {
pp = TMath::Sqrt(vect[2]*vect[2] + pt);
pt = TMath::Sqrt(pt);
ct = vect[2]/pp;
st = pt/pp;
cf = vect[0]/pt;
sf = vect[1]/pt;
}
rmx[0][0] = (cp * cf * ct) + (sp * sf);
rmx[0][1] = (cp * sf * ct) - (sp * cf);
rmx[0][2] = -(cp * st);
rmx[1][0] = -(cp * sf) + (sp * cf * ct);
rmx[1][1] = (cp * cf) + (sp * sf * ct);
rmx[1][2] = -(sp * st);
rmx[2][0] = (cf * st);
rmx[2][1] = (sf * st);
rmx[2][2] = ct;
}
//#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.