//////////////////////////////////////////////////////////////////////////
// //
// ATLFast MuonMaker class. //
// //
//////////////////////////////////////////////////////////////////////////
#include <TMCParticle.h>
#include <TRandom.h>
#include <TSystem.h>
#include <TFile.h>
#include <TDirectory.h>
#include "ATLFast.h"
#include "ATLFMCMaker.h"
#include "ATLFClusterMaker.h"
#include "ATLFMuonMaker.h"
#include "ATLFCluster.h"
#include "ATLFMuon.h"
const Float_t kPI = TMath::Pi();
const Float_t k2PI = 2*kPI;
const Int_t kMAXMUO = 100;
ClassImp(ATLFMuonMaker)
//_____________________________________________________________________________
ATLFMuonMaker::ATLFMuonMaker()
{
m_Nmuons = 0;
}
//_____________________________________________________________________________
ATLFMuonMaker::ATLFMuonMaker(const char *name, const char *title)
:ATLFMaker(name,title)
{
SetMinPT();
SetMaxEta();
SetRconeSep();
SetRconeDep();
SetMaxEdep();
m_Fruits = new TClonesArray("ATLFMuon",2, kFALSE);
m_BranchName = "Muons";
m_Nmuons = 0;
Save();
}
//_____________________________________________________________________________
ATLFMuonMaker::~ATLFMuonMaker()
{
//dummy
}
//_____________________________________________________________________________
void ATLFMuonMaker::AddMuon(Int_t code, Int_t mcparticle, Int_t mother, Int_t use, Int_t isol, Float_t eta, Float_t phi, Float_t pt, Int_t trigger)
{
// Add a new muon to the list of muons
//Note the use of the "new with placement" to create a new muon object.
//This complex "new" works in the following way:
// muons[i] is the value of the pointer for muon number i in the TClonesArray
// if it is zero, then a new muon must be generated. This typically
// will happen only at the first event
// If it is not zero, then the already existing object is overwritten
// by the new muon parameters.
// This technique should save a huge amount of time otherwise spent
// in the operators new and delete.
TClonesArray &muons = *(TClonesArray*)m_Fruits;
new(muons[m_Nmuons++]) ATLFMuon(code,mcparticle,mother,use,isol,eta,phi,pt,trigger);
}
//_____________________________________________________________________________
void ATLFMuonMaker::Clear(Option_t *option)
{
// Reset Muon Maker
m_Nmuons = 0;
ATLFMaker::Clear(option);
}
//_____________________________________________________________________________
void ATLFMuonMaker::Init()
{
// Create histograms and prepare resolution calculations - i.e load
// and prepare resolution histograms.
m_Mult = new TH1F("MuoMult","//muon multiplicity NOISOLATED",10,0,10);
m_MultIsol = new TH1F("MuoMultIsol","//muon multiplicity ISOLATED",10,0,10);
m_MultHard = new TH1F("MuoMultHard","muon multiplicity hard",10,0,10);
m_MultHardIsol = new TH1F("MuoMultHardIsol","muon multiplicity hard + isol",10,0,10);
m_PT = new TH1F("MuoPT","(pT- pTcru) vers pT muon ISOLATED",50,0,100);
m_Eta = new TH1F("MuoEta","(eta-etacru) vers pT muon ISOLATED",50,0,100);
m_Phi = new TH1F("MuoPhi","(phi-phicru) vers pT muon ISOLATED",50,0,100);
m_Counter = new TH1F("MuoCounter","counter vers pT muon ISOLATED",50,0,100);
m_Mass2mu = new TH1F("MuoMass2mu","mumu mass",50,80,100);
m_Mass2musubst = new TH1F("MuoMass2musubst","mumu mass subtract width",50,80,100);
m_Mass4mu = new TH1F("MuoMass4mu","4mu mass",50,120,140);
//translation of SUBROUTINE INI_RESOLUMU
// if "tracker only" option selected we don't need resolution histograms
if (gATLFast->SmearMuonOpt()==2) return;
// else: connect file with resolution histograms and process them...
Text_t *resfilename=0;
//to be portable between Unix and NT (slashes/backslahes problem)
//was strcpy(resfilename,gSystem->ExpandPathName("$ATLFAST/data/resolumu.root"));
resfilename = gSystem->ConcatFileName(gSystem->Getenv("ATLFAST"),"/data/resolumu.root");
if (gSystem->AccessPathName(resfilename,kReadPermission)){
Error("Init","Cannot open a file with Resolution histograms:n $ATLFAST/data/resolumu.root. n Please check whether this file/path exist");
gSystem->Exit(1,1);
}
m_Resms=new TH2F("m_Resms","",m_nPhi,0,m_nPhi,m_nEta,0,m_nEta);
m_Respr=new TH2F("m_Respr","",m_nPhi,0,m_nPhi,m_nEta,0,m_nEta);
m_Refms=new TH2F("m_Refms","",m_nPhi,0,m_nPhi,m_nEta,0,m_nEta);
m_Refpr=new TH2F("m_Refpr","",m_nPhi,0,m_nPhi,m_nEta,0,m_nEta);
// prevent these histograms from being added to a list of ATLFMuonMaker
// histograms - these are only temporary histograms used in calculations,
// and thus needn't be saved .
m_Resms->SetDirectory(0);
m_Respr->SetDirectory(0);
m_Refms->SetDirectory(0);
m_Refpr->SetDirectory(0);
// prepare File environment, because we'll open new file now...
TDirectory *dirsav=gDirectory;
// get resolution histograms form file:
m_file=new TFile(resfilename,"READ");
delete [] resfilename;
TH2F *rgresm=(TH2F*)m_file->Get("rgresm");
TH2F *rtresm=(TH2F*)m_file->Get("rtresm");
TH2F *rgrefm=(TH2F*)m_file->Get("rgrefm");
TH2F *rtrefm=(TH2F*)m_file->Get("rtrefm");
// close the file and restore old File environment:
m_file->Close();
dirsav->cd();
// hard code STD Pt and FEET Pt values - we won't keep it in file anymore
Float_t pt1s=100.0;
Float_t pt2s=1000.0;
Float_t pt1f=100.0;
Float_t pt2f=1000.0;
m_dEta = 2.7/(Float_t)(m_nEta-1);
m_minEta = -m_dEta/2.0;
m_maxEta = 2.7 + m_dEta/2.0;
m_dPhi = 45.0 / (Float_t)(m_nPhi-1);
m_minPhi = -m_dPhi/2.0;
m_maxPhi = 45.0 + m_dPhi/2.0;
Int_t leta=0;
Int_t lphi;
Int_t ieta;
Int_t iphi;
Float_t eta,delos1s,delos2s,difloss,delos1f,delos2f,diflosf;
Float_t rgres,rtres,rgref,rtref,amulsca,bmulsca;
for (ieta=1;ieta<=m_nEta;ieta++){
leta++;
eta = m_dEta*(Float_t)(ieta-1);
delos1s = Delosmu(pt1s,eta);
delos2s = Delosmu(pt2s,eta);
difloss = delos2s*delos2s - delos1s*delos1s;
delos1f = Delosmu(pt1f,eta);
delos2f = Delosmu(pt2f,eta);
diflosf = delos2f*delos2f - delos1f*delos1f;
// now recalculate calibration histograms - we store the results
// in our temporary histograms. These histograms are used in
// ResolMuo() method...
lphi=0;
for(iphi=1;iphi<=m_nPhi;iphi++){
lphi++;
rgres=rgresm->GetCellContent(lphi,leta)/100.0;
rtres=rtresm->GetCellContent(lphi,leta)/100.0;
rgref=rgrefm->GetCellContent(lphi,leta)/100.0;
rtref=rtrefm->GetCellContent(lphi,leta)/100.0;
if (rgres<1.0){
bmulsca = (rtres*rtres - rgres*rgres -difloss)/
( pt2s*pt2s - pt1s*pt1s);
amulsca = rgres*rgres-bmulsca*pt1s*pt1s-delos1s*delos1s;
} else {
amulsca = 1;
bmulsca = 0;
}
m_Resms->SetCellContent(lphi,leta,amulsca);
m_Respr->SetCellContent(lphi,leta,amulsca);
if (rgref<1.0){
bmulsca = (rtref*rtref - rgref*rgref -diflosf)/
(pt2f*pt2f - pt1f*pt1f);
amulsca = rgref*rgref-bmulsca*pt1f*pt1f-delos1f*delos1f;
} else {
amulsca = 1;
bmulsca = 0;
}
m_Refms->SetCellContent(lphi,leta,amulsca);
m_Refpr->SetCellContent(lphi,leta,bmulsca);
}
}
}
//_____________________________________________________________________________
void ATLFMuonMaker::Finish()
{
// Function called by ATLFast::Finish at the end of the job
// Normalize histograms
m_PT->Divide(m_Counter);
m_Eta->Divide(m_Counter);
m_Phi->Divide(m_Counter);
}
//_____________________________________________________________________________
Int_t ATLFMuonMaker::Make()
{
//.............................................
// This function searches for muons, by scanning through
// the list of MC particles. If a muon is found, its momentum is
// smeared using function RESMUO.
// Three options for muon-momentum smearing are available:
// stand-alone Muon System
// Inner Detector alone
// combined
// The parametrization for the momentum smearing is coded in RESMUO.
// Isolated and non-isolated muons are stored in The TClonesArray of muons
// and the energy clusters associated with them are removed.
// Muons outside the ETA-coverage or below the p_T-threshold are lost.
//.............................................
Int_t muCode[kMAXMUO], muPart[kMAXMUO], muFlag[kMAXMUO], muMother[kMAXMUO], muIsol[kMAXMUO];
Int_t muxCode[kMAXMUO], muxPart[kMAXMUO], muxFlag[kMAXMUO], muxMother[kMAXMUO], muxIsol[kMAXMUO];
Float_t muEta[kMAXMUO], muPhi[kMAXMUO], muPT[kMAXMUO];
Float_t muxEta[kMAXMUO], muxPhi[kMAXMUO], muxPT[kMAXMUO];
Int_t sortindex[kMAXMUO];
Int_t nmuons = 0;
Int_t nmuonsx = 0;
Int_t i, k, ic, ij, KF;
//........................................................
//... Get pointers to MCparticles arrays and ClonesArray
//........................................................
ATLFMCMaker *mcarlo = gATLFast->MCMaker();
TClonesArray *particles = mcarlo->Fruits();
TMCParticle *part, *parent;
Float_t *pPT = mcarlo->PT();
Float_t *pEta = mcarlo->Eta();
Float_t *pPhi = mcarlo->Phi();
Int_t *Index= mcarlo->Index();
//........................................................
//... Get pointers to clustermaker arrays and ClonesArray
//........................................................
ATLFClusterMaker *clustermaker = gATLFast->ClusterMaker();
TClonesArray *clusters = clustermaker->Fruits();
ATLFCluster *cluster;
Float_t *cellEta = clustermaker->CellEta();
Float_t *cellPhi = clustermaker->CellPhi();
Float_t *cellSumPT = clustermaker->CellSumPT();
Int_t nclusters = clustermaker->Nclusters();
Int_t ncells = clustermaker->Ncells();
//........................................................
//....look for isolated muons, sort clusters common
//........................................................
Bool_t isolated;
Float_t ptcru, etacru, phicru;
Float_t sigma, pxmuo, pymuo, pzmuo, eemuo;
Float_t pt, eta, phi, atheta, abseta, deta, dphi;
Float_t adphi, ddr, edep;
Float_t theta = 0;
Int_t nstop = mcarlo->Nstop();
Float_t RconeSep2 = m_RconeSep*m_RconeSep;
Float_t RconeDep2 = m_RconeDep*m_RconeDep;
for (Int_t ind=0;ind<mcarlo->Nstables();ind++) {
i = Index[ind];
if (i < mcarlo->Nstart()) continue;
part = (TMCParticle*)particles->UncheckedAt(i);
KF = part->GetKF();
if (TMath::Abs(KF) != 13) continue;
//-- analyse non-decayed muons
isolated = kTRUE;
pt = pPT[i];
eta = pEta[i];
phi = pPhi[i];
ptcru = pt;
etacru = eta;
phicru = phi;
//.... apply smearing
if (gATLFast->Smearing() == 1) {
sigma = MuonResolution(gATLFast->SmearMuonOpt(),gATLFast->SmearMuonFun(), pt, eta, phi);
pxmuo = part->GetPx()/(1+sigma);
pymuo = part->GetPy()/(1+sigma);
pzmuo = part->GetPz()/(1+sigma);
eemuo = part->GetEnergy()/(1+sigma);
if (TMath::Abs(pzmuo/eemuo) < 1) {
theta = TMath::ACos(pzmuo/eemuo);
} else {
if (pzmuo > 0) theta = 0;
if (pzmuo < 0) theta = kPI;
}
atheta = TMath::Abs(TMath::Tan(0.5*theta));
if (atheta < 0.0001) atheta = 0.0001;
eta = -TMath::Log(atheta);
phi = mcarlo->Angle(pxmuo,pymuo);
pt = TMath::Sqrt(pxmuo*pxmuo + pymuo*pymuo);
}
if (pt < m_MinPT) continue;
abseta = TMath::Abs(eta);
if (abseta > m_MaxEta) continue;
//........................................................
//....check if muon isolated from clusters
//........................................................
for (k=0;k<nclusters;k++) {
cluster = (ATLFCluster*)clusters->UncheckedAt(k);
deta = eta - cluster->Eta();
dphi = phi - cluster->Phi();
adphi = TMath::Abs(dphi) - k2PI;
if (TMath::Abs(dphi) > kPI) ddr = deta*deta +adphi*adphi;
else ddr = deta*deta + dphi*dphi;
if (ddr < RconeSep2) isolated = kFALSE;
}
//........................................................
//....check on energy deposition of cells EDEP in cone RDEP
//........................................................
edep = 0;
for (ic=0;ic<ncells;ic++) {
deta = eta - cellEta[ic];
dphi = phi - cellPhi[ic];
adphi = TMath::Abs(dphi) - k2PI;
if (TMath::Abs(dphi) > kPI) ddr = deta*deta + adphi*adphi;
else ddr = deta*deta + dphi*dphi;
if (ddr < RconeDep2) edep += cellSumPT[ic];
}
if (edep > m_MaxEdep) isolated = kFALSE;
parent = (TMCParticle*)particles->UncheckedAt(part->GetParent()-1);
if (isolated) { // fill isolated muon structure
muCode[nmuons] = KF;
muPart[nmuons] = i;
muMother[nmuons] = parent->GetKF();
muFlag[nmuons] = 0;
muIsol[nmuons] = 1;
muEta[nmuons] = eta;
muPhi[nmuons] = phi;
muPT[nmuons] = pt;
nmuons++;
m_PT->Fill(pt, TMath::Abs(pt-ptcru));
m_Eta->Fill(pt,TMath::Abs(eta-etacru));
m_Phi->Fill(pt,TMath::Abs(phi-phicru));
m_Counter->Fill(pt);
} else { // fill non-isolated muon structure
muxCode[nmuonsx] = KF;
muxPart[nmuonsx] = i;
muxMother[nmuonsx] = parent->GetKF();
muxFlag[nmuonsx] = 1;
muxIsol[nmuonsx] = 0;
muxEta[nmuonsx] = eta;
muxPhi[nmuonsx] = phi;
muxPT[nmuonsx] = pt;
nmuonsx++;
}
} // end of loop on primary particles
m_MultIsol->Fill(nmuons+0.01);
m_Mult->Fill(nmuonsx+0.01);
//........................................................
//.....arrange isolated muons in falling E_T sequence
//........................................................
Int_t muonCode[kMAXMUO], muonPart[kMAXMUO], muonFlag[kMAXMUO], muonMother[kMAXMUO], muonIsol[kMAXMUO];
Float_t muonEta[kMAXMUO], muonPhi[kMAXMUO], muonPT[kMAXMUO];
gATLFast->SortDown(nmuons, muPT, sortindex);
for (i=0;i<nmuons;i++) {
ij = sortindex[i];
muonCode[i] = muCode[ij];
muonPart[i] = muPart[ij];
muonMother[i] = muMother[ij];
muonFlag[i] = muFlag[ij];
muonIsol[i] = muIsol[ij];
muonEta[i] = muEta[ij] ;
muonPhi[i] = muPhi[ij];
muonPT[i] = muPT[ij];
}
//........................................................
//.....arrange non-isolated muons in falling E_T sequence
//........................................................
Int_t muonxCode[kMAXMUO], muonxPart[kMAXMUO], muonxFlag[kMAXMUO], muonxMother[kMAXMUO], muonxIsol[kMAXMUO];
Float_t muonxEta[kMAXMUO], muonxPhi[kMAXMUO], muonxPT[kMAXMUO];
gATLFast->SortDown(nmuonsx, muxPT, sortindex);
for (i=0;i<nmuonsx;i++) {
ij = sortindex[i];
muonxCode[i] = muxCode[ij];
muonxPart[i] = muxPart[ij];
muonxMother[i] = muxMother[ij];
muonxFlag[i] = muxFlag[ij];
muonxIsol[i] = muxIsol[ij];
muonxEta[i] = muxEta[ij] ;
muonxPhi[i] = muxPhi[ij];
muonxPT[i] = muxPT[ij];
}
Float_t pxmuo1, pymuo1, pzmuo1, eemuo1;
Float_t pxmuo2, pymuo2, pzmuo2, eemuo2;
Float_t px4mu, py4mu, pz4mu, ee4mu;
Float_t amumu, ammc, ammcc;
//........................................................
//......mumu mass distribution
//........................................................
if (nmuons == 2 && muonIsol[0] == 1 && muonIsol[1] == 1 ) {
pxmuo1 = muonPT[0]*TMath::Cos(muonPhi[0]);
pymuo1 = muonPT[0]*TMath::Sin(muonPhi[0]);
pzmuo1 = muonPT[0]*TMath::SinH(muonEta[0]);
eemuo1 = muonPT[0]*TMath::CosH(muonEta[0]);
pxmuo2 = muonPT[1]*TMath::Cos(muonPhi[1]);
pymuo2 = muonPT[1]*TMath::Sin(muonPhi[1]);
pzmuo2 = muonPT[1]*TMath::SinH(muonEta[1]);
eemuo2 = muonPT[1]*TMath::CosH(muonEta[1]);
amumu = (eemuo1+eemuo2)*(eemuo1+eemuo2) - (pxmuo1+pxmuo2)*(pxmuo1+pxmuo2)
-(pymuo1+pymuo2)*(pymuo1+pymuo2) - (pzmuo1+pzmuo2)*(pzmuo1+pzmuo2);
if (amumu > 0) amumu = TMath::Sqrt(amumu);
m_Mass2mu->Fill(amumu);
//.....substract natural width
for (k=0;k<nstop;k++) {
part = (TMCParticle*)particles->UncheckedAt(k);
parent = (TMCParticle*)particles->UncheckedAt(part->GetParent()-1);
KF = part->GetKF();
if (KF == 13 && parent->GetKF() == 23) {
ammc = 91.173;
ammcc = parent->GetMass();
m_Mass2musubst->Fill(amumu-(ammcc-ammc));
}
}
}
//........................................................
//........mass of 4mu....cross-check with partons
//........................................................
if (nmuons == 4 && muonIsol[0] == 1 && muonIsol[1] == 1 && muonIsol[2] == 1 && muonIsol[3] == 1) {
px4mu = py4mu = pz4mu = ee4mu = 0;
for (k=0;k<4;k++) {
px4mu += muonPT[k]*TMath::Cos(muonPhi[k]);
py4mu += muonPT[k]*TMath::Sin(muonPhi[k]);
pz4mu += muonPT[k]*TMath::SinH(muonEta[k]);
ee4mu += muonPT[k]*TMath::CosH(muonEta[k]);
}
amumu = ee4mu*ee4mu -px4mu*px4mu -py4mu*py4mu - pz4mu*pz4mu;
if (amumu > 0) amumu = TMath::Sqrt(amumu);
m_Mass4mu->Fill(amumu);
}
//........................................................
//....cross-check with partons
//........................................................
Int_t imuo = 0;
Int_t imuoiso = 0;
Float_t ener;
Int_t KF2;
TMCParticle *part2;
for (i=0;i<nstop;i++) {
part = (TMCParticle*)particles->UncheckedAt(i);
KF = part->GetKF();
if (TMath::Abs(KF) == 13) {
abseta = TMath::Abs(pEta[i]);
ener = 0;
isolated = kTRUE;
for (k=0;k<nstop;k++) {
part2 = (TMCParticle*)particles->UncheckedAt(k);
KF2 = TMath::Abs(part2->GetKF());
if ((KF2 <= 21) && (i != k) && (KF2 != 12) && (KF2 != 14) && (KF2 !=16)) {
deta = pEta[i] - pEta[k];
dphi = pPhi[i] - pPhi[k];
adphi = TMath::Abs(dphi) - k2PI;
if (TMath::Abs(dphi) > kPI) ddr = deta*deta + adphi*adphi;
else ddr = deta*deta + dphi*dphi;
if (ddr < RconeSep2 && pPT[k] > clustermaker->MinETCalo()) isolated = kFALSE;
if (ddr < RconeDep2 && pPT[k] < clustermaker->MinETCalo()) ener += pPT[k];
}
}
if (ener > m_MaxEdep) isolated = kFALSE;
if (abseta < m_MaxEta && pPT[i] > m_MinPT) {
imuo++;
if (isolated) imuoiso++;
}
}
}
m_MultHard->Fill(imuo+0.01);
m_MultHardIsol->Fill(imuoiso+0.01);
//........................................................
//....Store muons in Root ClonesArray
//........................................................
// AddMuon(Int_t code, Int_t mcparticle, Int_t mother, Int_t use, Int_t isol, Float_t eta, Float_t phi, Float_t pt, Int_t trigger)
m_Nmuons = 0;
for (k=0;k<nmuons;k++) {
AddMuon(muonCode[k],
muonPart[k],
muonMother[k],
muonFlag[k],
muonIsol[k],
muonEta[k],
muonPhi[k],
muonPT[k],
1);
}
for (k=0;k<nmuonsx;k++) {
AddMuon(muonxCode[k],
muonxPart[k],
muonxMother[k],
muonxFlag[k],
muonxIsol[k],
muonxEta[k],
muonxPhi[k],
muonxPT[k],
1);
}
return 0;
}
//_____________________________________________________________________________
Float_t ATLFMuonMaker::MuonResolution(Int_t keymuo, Int_t keyfun, Float_t pt, Float_t eta, Float_t phi)
{
// parametrizes muon resolution for
// keyfun=1 parametrization from C. Guyot
// keyfun=2 parametrization from M.Virchaux & L.Chevalier
// keymuo=1 -- muon system stand-alone
// keymuo=2 -- only tracker
// keymuo=3 -- combined
// Till now, only simpler version : keyfun=2 is supported...
Float_t sigma = 0;
Float_t aa, bb, atrackpt;
Float_t sigmu,sigmuon,atrak,wmuon,sigtrak,wtrak,sigtr,corr,wtot;
Float_t abseta = TMath::Abs(eta);
Float_t absetapow = TMath::Power(abseta,10)/7000.;
if(keymuo == 1) {
while(1) {
gRandom->Rannor(aa,bb);
sigma = 0.01*ResolMuo(keyfun,abseta,phi,pt);
sigma = aa*sigma;
if(1.+sigma > 0) return sigma;
}
} else if(keymuo == 2) {
while(1) {
gRandom->Rannor(aa,bb);
atrak = 0.0005;
atrak = atrak*(1.+absetapow);
sigma = atrak*pt*aa + 0.012*bb;
sigma = sigma*(1.+absetapow);
if(1.+sigma > 0) return sigma;
}
} else if(keymuo == 3) {
while(1) {
gRandom->Rannor(aa,bb);;
sigmuon = 0.01*ResolMuo(keyfun,abseta,phi,pt);
sigmu = aa*sigmuon;
if(1.+sigmu > 0) break;
}
while(1) {
gRandom->Rannor(aa,bb);
atrak = 0.0005;
atrak = atrak*(1.+absetapow);
atrackpt= atrak*pt;
sigtrak = TMath::Sqrt(atrackpt*atrackpt + 0.012*0.012);
sigtr = atrak*pt*aa + 0.012*bb;
if(1.+sigtr > 0) break;
}
wmuon = 1.0/(sigmuon*sigmuon);
wtrak = 1.0/(sigtrak*sigtrak);
wtot = wmuon+wtrak;
corr = (wmuon*(1.0+sigmu)+wtrak*(1.0+sigtr))/wtot;
sigma = corr-1.0;
}
return sigma;
}
//_____________________________________________________________________________
Float_t ATLFMuonMaker::ResolMuo(Int_t keyfun, Float_t &pt, Float_t &eta, Float_t &phi)
{
// This is the implementation of RESOLUMU subroutine of fortran version.
//
// SUBROUTINE RESOLUMU(PT,ETA,PHI,RESOL)
//*-----------------------------------------------------------------------
//*
//* Calculates the stand-alone muon spectrometer resolution from
//* an interpolation between tables at pt=100 and 1000 GeV read from
//* the file 'resolumu.root' .
//* Spectro Configuration as for the TDR
//*
//* input: Pt in GeV, phi in degrees, eta
//* output: resol in percent (if resol ge.100. : no measurement)
//*-----------------------------------------------------------------------
//
// WARNING: only keyfun=2 is valid in current version.
//
if (keyfun!=2){
Error("ResolMuo"," keyfun=%i is invalid ! n Only keyfun=2 is valid in this version of ATLFASTn",keyfun);
gSystem->Exit(-1,1);
}
Float_t respr2,resms2,resol2,deloss;
Int_t leta,lphi;
Float_t resol = 1.0;
Float_t aeta = TMath::Abs(eta);
if (aeta < m_maxEta)
{
Float_t phid = phi;
if (phid < 0) phid += 360.0;
Int_t isec = (Int_t)(phid/45.0 + 1.0);
Float_t phis = phid-(isec-1)*45.0;
leta = (Int_t)((aeta-m_minEta)/m_dEta) +1;
lphi = (Int_t)((phis-m_minPhi)/m_dPhi) +1;
if (isec!=6 && isec!=7){ // standard sector
respr2 = m_Respr->GetCellContent(lphi,leta)*pt*pt;
resms2 = m_Resms->GetCellContent(lphi,leta);
} else {
if (isec == 6) lphi = m_nPhi-lphi+1;
respr2 = m_Refpr->GetCellContent(lphi,leta)*pt*pt;
resms2 = m_Refms->GetCellContent(lphi,leta);
}
if (resms2 < 1.) {
resol2 = resms2+respr2;
deloss = Delosmu(pt,eta);
resol = TMath::Sqrt(resol2+deloss*deloss);
}
}
resol *= 100.0;
return resol;
}
//_____________________________________________________________________________
Float_t ATLFMuonMaker::Delosmu(Float_t &pt, Float_t &eta)
{
// Contribution from energy losses fluctuations in the calos.
// This function is used internaly by ResolMuo() method.
Float_t fact;
const Float_t kFactFwd=TMath::Sqrt(14.0/10.5);
Float_t aeta = TMath::Abs(eta);
Float_t theta= 2.0*TMath::ATan(TMath::Exp(-aeta));
Float_t sinth= TMath::Sin(theta);
Float_t emu = TMath::Abs(pt) /sinth;
if (aeta<1.1) fact = 1.0/TMath::Sqrt(sinth); // barrrel;
else fact = kFactFwd; // Forward region;
return (0.480 +0.105*emu/100.0)/emu*fact;
}
//_____________________________________________________________________________
void ATLFMuonMaker::PrintInfo()
{
ATLFMaker::PrintInfo();
}
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.