//////////////////////////////////////////////////////////////////////////
// //
// ATLFast JetMaker class. //
// //
//////////////////////////////////////////////////////////////////////////
#include <TMCParticle.h>
#include "ATLFJetMaker.h"
#include "ATLFJet.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 kMAXJETS = 100;
ClassImp(ATLFJetMaker)
//_____________________________________________________________________________
ATLFJetMaker::ATLFJetMaker()
{
m_Njets = 0;
}
//_____________________________________________________________________________
ATLFJetMaker::ATLFJetMaker(const char *name, const char *title)
:ATLFMaker(name,title)
{
// Default Setters for jets
SetMinET();
SetEtaCoverage();
// Default Setters for BJets
SetBquarkMinPT();
SetBquarkMaxEta();
SetBRconeMatch();
// Default Setters for CJets
SetCquarkMinPT();
SetCquarkMaxEta();
SetCRconeMatch();
// Default Setters for TauJets
SetTMinPTHad();
SetTMaxEta();
SetTRconeMatch();
SetTRatioPT();
m_Fruits = new TClonesArray("ATLFJet",15, kFALSE);
m_BranchName = "Jets";
m_Njets = 0;
Save();
}
//_____________________________________________________________________________
ATLFJetMaker::~ATLFJetMaker()
{
//dummy
}
//_____________________________________________________________________________
void ATLFJetMaker::AddJet(Int_t code, Int_t ncells, Int_t nparticles,Int_t part,
Float_t eta0, Float_t phi0, Float_t eta, Float_t phi, Float_t pt)
{
// Add a new jet to the list of jets
//Note the use of the "new with placement" to create a new jet object.
//This complex "new" works in the following way:
// jets[i] is the value of the pointer for jet number i in the TClonesArray
// if it is zero, then a new jet must be generated. This typically
// will happen only at the first events
// If it is not zero, then the already existing object is overwritten printf("muons-jestem tutajn");
// by the new jet parameters.
// This technique should save a huge amount of time otherwise spent
// in the operators new and delete.
TClonesArray &jets = *(TClonesArray*)m_Fruits;
new(jets[m_Njets++]) ATLFJet(code,ncells,nparticles,part,eta0,phi0,eta,phi,pt);
}
//_____________________________________________________________________________
void ATLFJetMaker::Clear(Option_t *option)
{
// Reset Jet Maker
m_Njets = 0;
ATLFMaker::Clear(option);
}
//_____________________________________________________________________________
void ATLFJetMaker::Init()
{
// Jets histograms
m_Mult = new TH1F("JetMult","jet multiplicity",10,0,10);
m_Mass2jj = new TH1F("JetMass2jj","jetjet mass",50,0,150);
m_ERes4050 = new TH1F("JetERes4050","energy resolution PT=40-50 GeV",50,-0.5,0.5);
m_ERes200250 = new TH1F("JetERes200250","energy resolution PT=200-250 GeV",50,-0.5,0.5);
// BJets histograms
m_BMult = new TH1F("JetBMult","b-jets multiplicity",10,0,10);
m_BMultHard = new TH1F("JetBMultHard","b-quarks HARD multiplicity",10,0,10);
m_BPTFSR = new TH1F("JetBPTFSR","p_T b-jet/b-quark FSR",50,0,1.5);
m_BPTHard = new TH1F("JetBPTHard","p_T b-jet/b-quark HARD",50,0,1.5);
m_BEtaFSR = new TH1F("JetBEtaFSR","eta b-jet/b-quark FSR",50,0.5,1.5);
m_BEtaHard = new TH1F("JetBEtaHard","eta b-jet/b-quark HARD",50,0.5,1.5);
m_BPhiFSR = new TH1F("JetBPhiFSR","phi b-jet/b-quark FSR",50,0.5,1.5);
m_BPhiHard = new TH1F("JetBPhiHard","phi b-jet/b-quark HARD",50,0.5,1.5);
m_BMjbjb = new TH1F("JetBMjbjb","b-jets Mjbjb",50,0,200);
m_BMbbFSR = new TH1F("JetBMbbFSR","b-quark Mbb (after FSR)",50,0,200);
m_BMbb = new TH1F("JetBMbb","M_bb b-jets/b-quark",50,0.5,1.5);
m_BRmin = new TH1F("JetBRmin","Rmin b-jet, b-quark HARD",50,0.0,20.0);
// CJets histograms
m_CMult = new TH1F("JetCMult","c-jets multiplicity",10,0,10);
m_CMultHard = new TH1F("JetCMultHard","c-quarks HARD multiplicity",10,0,10);
m_CPTFSR = new TH1F("JetCPTFSR","p_T c-jet/c-quark FSR",50,0,1.5);
m_CPTHard = new TH1F("JetCPTHard","p_T c-jet/c-quark HARD",50,0,1.5);
m_CEtaFSR = new TH1F("JetCEtaFSR","eta c-jet/c-quark FSR",50,0.5,1.5);
m_CEtaHard = new TH1F("JetCEtaHard","eta c-jet/c-quark HARD",50,0.5,1.5);
m_CPhiFSR = new TH1F("JetCPhiFSR","phi c-jet/c-quark FSR",50,0.5,1.5);
m_CPhiHard = new TH1F("JetCPhiHard","phi c-jet/c-quark HARD",50,0.5,1.5);
m_CMjcjc = new TH1F("JetCMjcjc","c-jets Mjcjc",50,0,200);
m_CMccFSR = new TH1F("JetCMccFSR","c-quark Mcc (after FSR)",50,0,200);
m_CMcc = new TH1F("JetCMcc","M_cc c-jets/c-quark",50,0.5,1.5);
m_CRmin = new TH1F("JetCRmin","Rmin c-jet, c-quark HARD",50,0.0,20.0);
// TauJets histograms
m_TauJetMult = new TH1F("JetTauJetMult","tau-jets multiplicity",10,0,10);
m_TausMult = new TH1F("JetTausMult","taus multiplicity",10,0,10);
}
//_____________________________________________________________________________
void ATLFJetMaker::Finish()
{
// Function called by ATLFast::Finish at the end of the job
}
//_____________________________________________________________________________
Int_t ATLFJetMaker::Make()
{
//.............................................
// This function searches for jets, by scanning through
// the list of Clusters. Each cluster momentum is smeared using function RESHAD.
// Jets outside the ETA-coverage or below the p_T-threshold are lost.
// Jets are labelled as b-jets, c-jets and tau-jets.
//.............................................
MakeJets();
MakeBJets();
MakeCJets();
MakeTauJets();
return 0;
}
//_____________________________________________________________________________
void ATLFJetMaker::MakeJets()
{
//.............................................
// This function searches for jets, by scanning through
// the list of clusters. Each cluster momentum is smeared using function
// Jets outside the ETA-coverage or below the p_T-threshold are lost.
//.............................................
Int_t jetCode[kMAXJETS], jetCells[kMAXJETS], jetPart[kMAXJETS];
Float_t jetEta0[kMAXJETS], jetPhi0[kMAXJETS];
Float_t jetEta[kMAXJETS], jetPhi[kMAXJETS], jetPT[kMAXJETS];
//........................................................
//... Get pointers to clustermaker arrays and ClonesArray
//........................................................
ATLFMCMaker *mcarlo = gATLFast->MCMaker();
ATLFClusterMaker *clustermaker = gATLFast->ClusterMaker();
TClonesArray *clusters = clustermaker->Fruits();
ATLFCluster *cluster;
Int_t nclusters = clustermaker->Nclusters();
Int_t i,k,ij;
Float_t rconeb, rconef;
rconeb = clustermaker->DeltaRBarrel1();
rconef = clustermaker->DeltaRForward1();
if (gATLFast->Mode() == 2) {
rconeb= clustermaker->DeltaRBarrel2();
rconef= clustermaker->DeltaRForward2();
}
//........................................................
//.....smear clusters energy
//........................................................
Float_t eeclu, pxclu, pyclu, pzclu;
Float_t etaclu, ptclu, phiclu, sigma, rcon, atheta;
Float_t eta, phi, et;
Float_t theta = 0;
if (gATLFast->Smearing() == 1) {
for (i=0;i<nclusters;i++) {
cluster = (ATLFCluster*)clusters->UncheckedAt(i);
eta = cluster->Eta();
phi = cluster->Phi();
et = cluster->ET();
eeclu = et*TMath::CosH(eta);
if (TMath::Abs(eta) <= clustermaker->BarrelForwardEta()) rcon = rconeb;
else rcon = rconef;
sigma = clustermaker->HadronResolution(gATLFast->Luminosity(),
eeclu,
eta,
et,
rcon);
pxclu = et*TMath::Cos(phi)*(1+sigma);
pyclu = et*TMath::Sin(phi)*(1+sigma);
pzclu = et*TMath::SinH(eta)*(1+sigma);
eeclu = et*TMath::CosH(eta)*(1+sigma);
if(TMath::Abs(pzclu/eeclu) < 1) theta = TMath::ACos(pzclu/eeclu);
else {
if (pzclu > 0) theta = 0;
if (pzclu < 0) theta = kPI;
}
atheta = TMath::Abs(TMath::Tan(0.5*theta));
if (atheta < 0.0001) atheta = 0.0001;
etaclu = -TMath::Log(atheta);
phiclu = mcarlo->Angle(pxclu,pyclu);
ptclu = TMath::Sqrt(pxclu*pxclu + pyclu*pyclu);
cluster->SetET(ptclu);
cluster->SetEta(etaclu);
cluster->SetPhi(phiclu);
if(ptclu > 40 && ptclu < 50) m_ERes4050->Fill(sigma/(1+sigma));
if(ptclu > 200 && ptclu < 250) m_ERes200250->Fill(sigma/(1+sigma));
}
}
//........................................................
//.....add nonisolated muons to jets
//........................................................
ATLFMuonMaker *muonmaker = gATLFast->MuonMaker();
TClonesArray *muons = muonmaker->Fruits();
ATLFMuon *muon;
Int_t nmuons = muonmaker->Nmuons();
Float_t adphi, ddr, deta, dphi, pt, dr;
Float_t cluEta, cluPhi, cluET;
Int_t muclu;
muclu=-1;
for (i=0;i<nmuons;i++) {
muon = (ATLFMuon*)muons->UncheckedAt(i);
if ( muon->Isolated() == 1) continue;
eta = muon->Eta();
phi = muon->Phi();
pt = muon->PT();
dr = 100;
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 < dr) {
muclu = k;
dr = ddr;
}
}
if(muclu <0) continue;
cluster = (ATLFCluster*)clusters->UncheckedAt(muclu);
cluEta = cluster->Eta();
cluPhi = cluster->Phi();
cluET = cluster->ET();
if ((TMath::Abs(cluEta) <= clustermaker->BarrelForwardEta() && dr < rconeb) ||
(TMath::Abs(cluEta) > clustermaker->BarrelForwardEta() && dr < rconef)) {
pxclu = cluET*TMath::Cos(cluPhi) + pt*TMath::Cos(phi);
pyclu = cluET*TMath::Sin(cluPhi) + pt*TMath::Sin(phi);
pzclu = cluET*TMath::SinH(cluEta) + pt*TMath::SinH(eta);
eeclu = cluET*TMath::CosH(cluEta) + pt*TMath::CosH(eta);
ptclu = TMath::Sqrt(pxclu*pxclu + pyclu*pyclu);
etaclu = mcarlo->Rapidity(ptclu, pzclu);
phiclu = mcarlo->Angle(pxclu, pyclu);
cluster->SetEta0(etaclu);
cluster->SetPhi0(phiclu);
cluster->SetEta(etaclu);
cluster->SetPhi(phiclu);
cluster->SetET(ptclu);
muon->SetUseFlag(0);
}
}
//........................................................
//.... store accepted jets in local arrays jetXXX and flag clusters
//........................................................
Int_t njets = 0;
for (k=0;k<nclusters;k++) {
cluster = (ATLFCluster*)clusters->UncheckedAt(k);
if (cluster->ET() > m_MinET && TMath::Abs(cluster->Eta()) < m_EtaCoverage) {
jetCode[njets] = cluster->KFcode();
jetCells[njets] = cluster->Ncells();
jetPart[njets] = cluster->Nparticles();
jetEta0[njets] = cluster->Eta0();
jetPhi0[njets] = cluster->Phi0();
jetEta[njets] = cluster->Eta();
jetPhi[njets] = cluster->Phi();
jetPT[njets] = cluster->ET();
njets++;
cluster->SetUseFlag(0);
}
}
m_Mult->Fill(njets+0.01);
//........................................................
//.....arrange jets in falling E_T sequence
//........................................................
Int_t sortindex[kMAXJETS];
Int_t jetsCode[kMAXJETS], jetsCells[kMAXJETS], jetsPart[kMAXJETS];
Float_t jetsEta0[kMAXJETS], jetsPhi0[kMAXJETS];
Float_t jetsEta[kMAXJETS], jetsPhi[kMAXJETS], jetsPT[kMAXJETS];
gATLFast->SortDown(njets, jetPT, sortindex);
for (i=0;i<njets;i++) {
ij = sortindex[i];
jetsCode[i] = jetCode[ij];
jetsCells[i] = jetCells[ij];
jetsPart[i] = jetPart[ij];
jetsEta0[i] = jetEta0[ij];
jetsPhi0[i] = jetPhi0[ij];
jetsEta[i] = jetEta[ij];
jetsPhi[i] = jetPhi[ij];
jetsPT[i] = jetPT[ij];
}
//........................................................
//.....jet-jet mass distribution
//........................................................
Float_t pxjet1, pyjet1, pzjet1, eejet1;
Float_t pxjet2, pyjet2, pzjet2, eejet2;
Float_t amjj;
if (njets == 2) {
pxjet1 = jetsPT[0]*TMath::Cos(jetsPhi[0]);
pyjet1 = jetsPT[0]*TMath::Sin(jetsPhi[0]);
pzjet1 = jetsPT[0]*TMath::SinH(jetsEta[0]);
eejet1 = jetsPT[0]*TMath::CosH(jetsEta[0]);
pxjet2 = jetsPT[1]*TMath::Cos(jetsPhi[1]);
pyjet2 = jetsPT[1]*TMath::Sin(jetsPhi[1]);
pzjet2 = jetsPT[1]*TMath::SinH(jetsEta[1]);
eejet2 = jetsPT[1]*TMath::CosH(jetsEta[1]);
amjj = (eejet1+eejet2)*(eejet1+eejet2) - (pxjet1+pxjet2)*(pxjet1+pxjet2)
-(pyjet1+pyjet2)*(pyjet1+pyjet2) - (pzjet1+pzjet2)*(pzjet1+pzjet2);
if (amjj> 0) amjj= TMath::Sqrt(amjj);
m_Mass2jj->Fill(amjj);
}
//........................................................
//....Store jets in Root ClonesArray
//........................................................
//AddJet(Int_t code, Int_t ncells, Int_t nparticles,Int_t part,
// Float_t etainit, Float_t phiinit, Float_t eta, Float_t phi, Float_t pt)
m_Njets = 0;
for (k=0;k<njets;k++) {
AddJet(jetsCode[k],
jetsCells[k],
jetsPart[k],
0,
jetsEta0[k],
jetsPhi0[k],
jetsEta[k],
jetsPhi[k],
jetsPT[k]);
}
}
//_____________________________________________________________________________
void ATLFJetMaker::MakeBJets()
{
//.............................................
// This function selects jets which can be associated to b-quarks,
// by scanning through the list of particles b-quarks(after FSR),
// which pass the acceptance cuts in PT and ETA are identified.
// For each of these b-quarks, a reconstructed jet associated with it
// is searched for. These jets are labelled as b-jets if certain acceptance
// criteria are fulfilled.
//.............................................
Int_t i, k, KF, nparticles;
//........................................................
//... Get pointers to MCparticles arrays and ClonesArray
//........................................................
ATLFMCMaker *mcarlo = gATLFast->MCMaker();
TClonesArray *particles = mcarlo->Fruits();
TMCParticle *part, *child;
nparticles = particles->GetEntriesFast();
Float_t *pPT = mcarlo->PT();
Float_t *pEta = mcarlo->Eta();
Float_t *pPhi = mcarlo->Phi();
//........................................................
//... Get pointers to Jets arrays and ClonesArray
//........................................................
ATLFJetMaker *jetmaker = gATLFast->JetMaker();
TClonesArray *jets = jetmaker->Fruits();
ATLFJet *jet;
Int_t njets = jets->GetEntriesFast();
Int_t njetb = 0;
Int_t jetb;
Bool_t bjet;
Float_t pz, pt, eta, phi, deta, dphi;
Float_t adphi, dr, ddr;
//........................................................
//....look for b-jets
//........................................................
for (i=mcarlo->Nstart();i<nparticles;i++) {
part = (TMCParticle*)particles->UncheckedAt(i);
KF = part->GetKF();
if (TMath::Abs(KF) != 5) continue;
//........................................................
//....if there is a b-quark found before hadronization
//....if there are still jets
//........................................................
if (njets <= 0) continue;
bjet = kTRUE;
jetb = -1;
//........................................................
//....and this b-quark is the last one in the FSR cascade
//........................................................
if (part->GetFirstChild()) {
for (k= part->GetFirstChild(); k<=part->GetLastChild();k++) {
child = (TMCParticle*)particles->UncheckedAt(k-1);
if (child && TMath::Abs(child->GetKF()) == 5) bjet = kFALSE;
}
}
// it is not filled, I don't know why???
// pt = pPT[i];
// eta = pEta[i];
// phi = pPhi[i];
pt=TMath::Sqrt(part->GetPx() * part->GetPx() + part->GetPy() * part->GetPy());
eta= mcarlo->Rapidity(pt,part->GetPz());
phi= mcarlo->Angle(part->GetPx(),part->GetPy());
if (pt < m_BquarkMinPT) bjet = kFALSE;
if (TMath::Abs(eta) > m_BquarkMaxEta) bjet = kFALSE;
if(!bjet)continue;
//........................................................
// === mark b-jet
//........................................................
dr = 100;
for (k=0;k<njets;k++) {
jet = (ATLFJet*)jets->UncheckedAt(k);
deta = eta - jet->Eta();
dphi = phi - jet->Phi();
adphi = TMath::Abs(dphi) - k2PI;
if (TMath::Abs(dphi) > kPI) ddr = TMath::Sqrt(deta*deta + adphi*adphi);
else ddr = TMath::Sqrt(deta*deta + dphi*dphi);
if (ddr < dr) jetb = k;
dr = TMath::Min(ddr, dr);
}
if (dr > m_BRconeMatch) {
bjet = kFALSE;
jetb = 0;
}
//........................................................
// === labell b-jet
//........................................................
if (bjet) {
jet = (ATLFJet*)jets->UncheckedAt(jetb);
jet->SetKFcode(5);
jet->SetPart(i);
m_BPTFSR->Fill(jet->PT()/pt);
if (eta !=0 ) m_BEtaFSR->Fill(jet->Eta()/eta);
if (phi != 0) m_BPhiFSR->Fill(jet->Phi()/phi);
njetb++;
}
}
m_BMult->Fill(njetb+0.01);
//........................................................
// === check mass resolution
//........................................................
Float_t pjx1,pjy1,pjz1,pje1;
Float_t pjx2,pjy2,pjz2,pje2;
Float_t E, px, py, massbb,massjbjb;
ATLFJet *jet1, *jet2;
TMCParticle *part1, *part2;
if (njetb == 2) {
for (i=0;i<njets-1;i++) {
jet1 = (ATLFJet*)jets->UncheckedAt(i);
if (TMath::Abs(jet1->KFcode()) != 5) continue;
pjx1 = jet1->PT()*TMath::Cos(jet1->Phi());
pjy1 = jet1->PT()*TMath::Sin(jet1->Phi());
pjz1 = jet1->PT()*TMath::SinH(jet1->Eta());
pje1 = jet1->PT()*TMath::CosH(jet1->Eta());
for (k=i+1;k<njets;k++) {
jet2 = (ATLFJet*)jets->UncheckedAt(k);
if (TMath::Abs(jet2->KFcode()) != 5) continue;
pjx2 = jet2->PT()*TMath::Cos(jet2->Phi());
pjy2 = jet2->PT()*TMath::Sin(jet2->Phi());
pjz2 = jet2->PT()*TMath::SinH(jet2->Eta());
pje2 = jet2->PT()*TMath::CosH(jet2->Eta());
massjbjb = (pje1+pje2)*(pje1+pje2) -(pjx1+pjx2)*(pjx1+pjx2)
-(pjy1+pjy2)*(pjy1+pjy2) -(pjz1+pjz2)*(pjz1+pjz2);
if (massjbjb > 0) massjbjb = TMath::Sqrt(massjbjb);
else massjbjb = 0;
m_BMjbjb->Fill(massjbjb);
part1 = (TMCParticle*)particles->UncheckedAt(jet1->Part());
part2 = (TMCParticle*)particles->UncheckedAt(jet2->Part());
E = part1->GetEnergy() + part2->GetEnergy();
px = part1->GetPx() + part2->GetPx();
py = part1->GetPy() + part2->GetPy();
pz = part1->GetPz() + part2->GetPz();
massbb = E*E -px*px -py*py -pz*pz;
if (massbb > 0) massbb = TMath::Sqrt(massbb);
else massbb = 0;
m_BMbbFSR->Fill(massbb);
if (massbb != 0) m_BMbb->Fill(massjbjb/massbb);
}
}
}
//........................................................
// === check partons
//........................................................
Int_t iquab = 0;
Int_t ibjet = -1;
Int_t nstop = mcarlo->Nstop();
for (i=0;i<nstop;i++) {
part = (TMCParticle*)particles->UncheckedAt(i);
KF = part->GetKF();
if (TMath::Abs(KF) != 5) continue;
if (TMath::Abs(pEta[i]) >= m_BquarkMaxEta || pPT[i] <= m_MinET) continue;
iquab++;
dr = 100;
for (k=0;k<njets;k++) {
jet = (ATLFJet*)jets->UncheckedAt(k);
if (TMath::Abs(jet->KFcode()) != 5) continue;
deta = pEta[i] - jet->Eta();
dphi = pPhi[i] - jet->Phi();
adphi = TMath::Abs(dphi) - k2PI;
if (TMath::Abs(dphi) > kPI) ddr = deta*deta + adphi*adphi;
else ddr = deta*deta + dphi*dphi;
if (ddr < dr) ibjet = k;
if (ddr < dr) dr = ddr;
}
if (ibjet < 0) continue; //<===========WARNING must be protected
jet = (ATLFJet*)jets->UncheckedAt(ibjet);
m_BRmin->Fill(TMath::Sqrt(dr));
m_BPTHard->Fill(jet->PT()/pPT[i]);
if (pEta[i] != 0) m_BEtaHard->Fill(jet->Eta()/pEta[i]);
if (pPhi[i] != 0) m_BPhiHard->Fill(jet->Phi()/pPhi[i]);
}
m_BMultHard->Fill(iquab+0.01);
}
//_____________________________________________________________________________
void ATLFJetMaker::MakeCJets()
{
//.............................................
// This function selects jets which can be associated to c-quarks,
// by scanning through the list of particles c-quarks(after FSR),
// which pass the acceptance cuts in PT and ETA are identified.
// For each of these c-quarks, a reconstructed jet associated with it
// is searched for. These jets are labelled as c-jets if certain acceptance
// criteria are fulfilled.
//.............................................
Int_t i, k, KF, nparticles;
//........................................................
//... Get pointers to MCparticles arrays and ClonesArray
//........................................................
ATLFMCMaker *mcarlo = gATLFast->MCMaker();
TClonesArray *particles = mcarlo->Fruits();
TMCParticle *part, *child;
nparticles = particles->GetEntriesFast();
Float_t *pPT = mcarlo->PT();
Float_t *pEta = mcarlo->Eta();
Float_t *pPhi = mcarlo->Phi();
//........................................................
//... Get pointers to Jets arrays and ClonesArray
//........................................................
ATLFJetMaker *jetmaker = gATLFast->JetMaker();
TClonesArray *jets = jetmaker->Fruits();
ATLFJet *jet;
Int_t njets = jets->GetEntriesFast();
Int_t njetc = 0;
Int_t jetc;
Bool_t cjet;
Float_t pz, pt, eta, phi, deta, dphi;
Float_t adphi, dr, ddr;
//........................................................
//....look for c-jets
//........................................................
for (i=mcarlo->Nstart();i<nparticles;i++) {
part = (TMCParticle*)particles->UncheckedAt(i);
KF = part->GetKF();
if (TMath::Abs(KF) != 4) continue;
//........................................................
//....if there is a c-quark found before hadronization
//....if there are still jets
//........................................................
if (njets <= 0) continue;
cjet = kTRUE;
jetc = -1;
//........................................................
//....and this c-quark is the last one in the FSR cascade
//........................................................
if (part->GetFirstChild()) {
for (k= part->GetFirstChild(); k<=part->GetLastChild();k++) {
child = (TMCParticle*)particles->UncheckedAt(k-1);
if (child && TMath::Abs(child->GetKF()) == 4) cjet = kFALSE;
}
}
// it is not filled, I don't know why???
// pt = pPT[i];
// eta = pEta[i];
// phi = pPhi[i];
pt=TMath::Sqrt(part->GetPx() * part->GetPx() + part->GetPy() * part->GetPy());
eta= mcarlo->Rapidity(pt,part->GetPz());
phi= mcarlo->Angle(part->GetPx(),part->GetPy());
if (pt < m_CquarkMinPT) cjet = kFALSE;
if (TMath::Abs(eta) > m_CquarkMaxEta) cjet = kFALSE;
//........................................................
// === mark c-jet
//........................................................
dr = 100;
for (k=0;k<njets;k++) {
jet = (ATLFJet*)jets->UncheckedAt(k);
deta = eta - jet->Eta();
dphi = phi - jet->Phi();
adphi = TMath::Abs(dphi) - k2PI;
if (TMath::Abs(dphi) > kPI) ddr = TMath::Sqrt(deta*deta + adphi*adphi);
else ddr = TMath::Sqrt(deta*deta + dphi*dphi);
if (ddr < dr) jetc = k;
dr = TMath::Min(ddr, dr);
}
if (dr > m_CRconeMatch) {
cjet = kFALSE;
jetc = 0;
}
//........................................................
// === labell c-jet
//........................................................
if (cjet) {
jet = (ATLFJet*)jets->UncheckedAt(jetc);
jet->SetKFcode(4);
jet->SetPart(i);
m_CPTFSR->Fill(jet->PT()/pt);
if (eta !=0 ) m_CEtaFSR->Fill(jet->Eta()/eta);
if (phi != 0) m_CPhiFSR->Fill(jet->Phi()/phi);
njetc++;
}
}
m_CMult->Fill(njetc+0.01);
//........................................................
// === check mass resolution
//........................................................
Float_t pjx1,pjy1,pjz1,pje1;
Float_t pjx2,pjy2,pjz2,pje2;
Float_t E, px, py, masscc,massjcjc;
ATLFJet *jet1, *jet2;
TMCParticle *part1, *part2;
if (njetc == 2) {
for (i=0;i<njets-1;i++) {
jet1 = (ATLFJet*)jets->UncheckedAt(i);
if (TMath::Abs(jet1->KFcode()) != 4) continue;
pjx1 = jet1->PT()*TMath::Cos(jet1->Phi());
pjy1 = jet1->PT()*TMath::Sin(jet1->Phi());
pjz1 = jet1->PT()*TMath::SinH(jet1->Eta());
pje1 = jet1->PT()*TMath::CosH(jet1->Eta());
for (k=i+1;k<njets;k++) {
jet2 = (ATLFJet*)jets->UncheckedAt(k);
if (TMath::Abs(jet2->KFcode()) != 4) continue;
pjx2 = jet2->PT()*TMath::Cos(jet2->Phi());
pjy2 = jet2->PT()*TMath::Sin(jet2->Phi());
pjz2 = jet2->PT()*TMath::SinH(jet2->Eta());
pje2 = jet2->PT()*TMath::CosH(jet2->Eta());
massjcjc = (pje1+pje2)*(pje1+pje2) -(pjx1+pjx2)*(pjx1+pjx2)
-(pjy1+pjy2)*(pjy1+pjy2) -(pjz1+pjz2)*(pjz1+pjz2);
if (massjcjc > 0) massjcjc = TMath::Sqrt(massjcjc);
else massjcjc = 0;
m_CMjcjc->Fill(massjcjc);
part1 = (TMCParticle*)particles->UncheckedAt(jet1->Part());
part2 = (TMCParticle*)particles->UncheckedAt(jet2->Part());
E = part1->GetEnergy() + part2->GetEnergy();
px = part1->GetPx() + part2->GetPx();
py = part1->GetPy() + part2->GetPy();
pz = part1->GetPz() + part2->GetPz();
masscc = E*E -px*px -py*py -pz*pz;
if (masscc > 0) masscc = TMath::Sqrt(masscc);
else masscc = 0;
m_CMccFSR->Fill(masscc);
if (masscc != 0) m_CMcc->Fill(massjcjc/masscc);
}
}
}
//........................................................
// === check partons
//........................................................
Int_t iquac = 0;
Int_t icjet = -1;
Int_t nstop = mcarlo->Nstop();
for (i=0;i<nstop;i++) {
part = (TMCParticle*)particles->UncheckedAt(i);
KF = part->GetKF();
if (TMath::Abs(KF) != 4) continue;
if (TMath::Abs(pEta[i]) >= m_CquarkMaxEta || pPT[i] <= m_MinET) continue;
iquac++;
dr = 18;
for (k=0;k<njets;k++) {
jet = (ATLFJet*)jets->UncheckedAt(k);
if (TMath::Abs(jet->KFcode()) != 4) continue;
deta = pEta[i] - jet->Eta();
dphi = pPhi[i] - jet->Phi();
adphi = TMath::Abs(dphi) - k2PI;
if (TMath::Abs(dphi) > kPI) ddr = deta*deta + adphi*adphi;
else ddr = deta*deta + dphi*dphi;
if (ddr < dr) icjet = k;
if (ddr < dr) dr = ddr;
}
if (icjet < 0) continue; //<===========WARNING must be protected
jet = (ATLFJet*)jets->UncheckedAt(icjet);
m_CRmin->Fill(TMath::Sqrt(dr));
m_CPTHard->Fill(jet->PT()/pPT[i]);
if (pEta[i] != 0) m_CEtaHard->Fill(jet->Eta()/pEta[i]);
if (pPhi[i] != 0) m_CPhiHard->Fill(jet->Phi()/pPhi[i]);
}
m_CMultHard->Fill(iquac+0.01);
}
//_____________________________________________________________________________
void ATLFJetMaker::MakeTauJets()
{
//.............................................
// This function selects jets which can be associated with TAU leptons
// decaying hadronically. The list of particles is scanned for hadronic
// decays of TAU leptons and if they pass acceptance criteria and
// associated jet is found, the jet is market as a TauJet.
//.............................................
Int_t i, k, ij, KF, nparticles, in;
//........................................................
//... Get pointers to MCparticles arrays and ClonesArray
//........................................................
ATLFMCMaker *mcarlo = gATLFast->MCMaker();
TClonesArray *particles = mcarlo->Fruits();
TMCParticle *part, *partin, *child;
nparticles = particles->GetEntriesFast();
//........................................................
//... Get pointers to Jets arrays and ClonesArray
//........................................................
ATLFJetMaker *jetmaker = gATLFast->JetMaker();
TClonesArray *jets = jetmaker->Fruits();
ATLFJet *jet;
Int_t njets = jets->GetEntriesFast();
Int_t njettau = 0;
Int_t ntau = 0;
Int_t jettau;
Bool_t tau;
Bool_t taujet;
Float_t pt, eta, phi, deta, dphi;
Float_t adphi, dr, ddr;
Float_t dpx, dpy, dpz;
//........................................................
//....look for tau-jets
//........................................................
for (i=mcarlo->Nstart();i<nparticles;i++) {
part = (TMCParticle*)particles->UncheckedAt(i);
KF = part->GetKF();
if (TMath::Abs(KF) != 15) continue;
//........................................................
//....if there is a tau-quark found before hadronization
//....if there are still jets
//........................................................
if (njets <= 0) continue;
taujet = kTRUE;
tau = kTRUE;
jettau = -1;
//........................................................
//....and this tau decay hadronically (tau-nu is in the cascade)
//........................................................
in = -1; // <===========Initialize
for (k= part->GetFirstChild(); k<=part->GetLastChild();k++) {
child = (TMCParticle*)particles->UncheckedAt(k-1);
if (!child) continue;
ij = child->GetKF();
if (TMath::Abs(ij) == 11 || TMath::Abs(ij) == 13) tau = kFALSE;
if (ij == 16) in = k-1;
}
if (in < 0) continue;
partin = (TMCParticle*)particles->UncheckedAt(in);
if (!partin) continue;
dpx = part->GetPx() -partin->GetPx();
dpy = part->GetPy() -partin->GetPy();
dpz = part->GetPz() -partin->GetPz();
pt = TMath::Sqrt(dpx*dpx + dpy*dpy);
eta = mcarlo->Rapidity(pt, dpz);
phi = mcarlo->Angle(dpx, dpy);
if (pt < m_TMinPTHad) tau= kFALSE;
if (TMath::Abs(eta) > m_TMaxEta) tau= kFALSE;
if (tau) ntau++;
if (!tau) continue;
//........................................................
// === mark tau-jet
//........................................................
dr = 100;
for (k=0;k<njets;k++) {
jet = (ATLFJet*)jets->UncheckedAt(k);
deta = eta - jet->Eta();
dphi = phi - jet->Phi();
adphi = TMath::Abs(dphi) - k2PI;
if (TMath::Abs(dphi) > kPI) ddr = TMath::Sqrt(deta*deta + adphi*adphi);
else ddr = TMath::Sqrt(deta*deta + dphi*dphi);
if (ddr < dr) jettau = k;
dr = TMath::Min(ddr, dr);
}
if (dr > m_TRconeMatch) {
taujet = kFALSE;
jettau = -1;
}
if (taujet) {
jet = (ATLFJet*)jets->UncheckedAt(jettau);
if (TMath::Abs(jet->Eta()) < 2.5 && pt/jet->PT() > m_TRatioPT) {
jet->SetKFcode(KF);
jet->SetPart(i);
njettau++;
}
}
}
m_TauJetMult->Fill(njettau+0.01);
m_TausMult->Fill(ntau+0.01);
}
//_____________________________________________________________________________
void ATLFJetMaker::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.