How to calculate missing four momentum ?

[Li Huang]

Hi, I am working on some e+ e- collider so I want to calculate the missing four momentum , not just the missing transverse energy which delphes result provide.
If I want to do like this: missing-momentum = init-momentum - detected-momentum.
so , how to get all the detected-momentum? shall I loop all the particle, jet , photon, etc...
Or, maybe you have some better way to get the missing four momentum.Can you tell me ?

[Michele Selvaggi]

Hi ,

you can use the ILD card (cards/delphes_card_ILD.tcl) which produces the MissingET branch by default.
You don't have to recalculate it. From there you can access the 3 components of the Missing Energy (MET,Eta,Phi) or (Px,Py,Pz).
You can access them in this way:

if(branchMissingET->GetEntriesFast() > 0)
{
   met = (MissingET*) branchMissingET->At(0);
   TVector3 V = met->P4().Vect();
   cout<<V.Mag()<<","<<V.Theta()<<","<<V.Phi()<<endl;  // E,Theta,Phi
   cout<<V.Px()<<","<<V.Py()<<","<<V.Pz()<<endl; //Px, Py, Pz
}

[Li Huang]

Hi,

Thanks for your answer, but I have a question,
Is the MET means missing transverse energy? I want, not only the missing transverse momentum but also the z direction momentum.

[Michele Selvaggi]

Hi,

We just called it MissingET since Delphes was initially developed for p-p collider where the "z" component has no meaning.
However, as shown in the small example above, the MissingET collection gives you also access to the "z" component (although this is something you should access only in the case of e+e- collision).

You have access to the 3x,y,z component as shown above

if(branchMissingET->GetEntriesFast() > 0)
{
   met = (MissingET*) branchMissingET->At(0);
   cout<<V.Px()<<","<<V.Py()<<","<<V.Pz()<<endl; //Px, Py, Pz
}

Cheers,
Michele

[Li Huang]

Hi,

Get it, Thanks a lot for your help!

[Li Huang]

Hi,

I try to generate a process e+ e- > z h , then go through pythia6 and Delphes with ILD card. Since there are 20% for z decay into neutrino , I guess I can use the missing four momentum to get the invariant mass of z. Then, I use,

if(branchMissingET->GetEntriesFast() > 0)
{

met = (MissingET*) branchMissingET->At(0);
TVector3 V = met->P4().Vect();
cout<<V.Mag()<<","<<V.Theta()<<","<<V.Phi()<<endl; E,Theta,Phi
cout<<V.Px()<<","<<V.Py()<<","<<V.Pz()<<endl; Px, Py, Pz
momentummet.SetPxPyPzE(V.Px(), V.Py(), V.Pz(),V.Mag());

}

Then I plot momentummet.M() , the invariant mass , but the results are almost 0 . There must be some thing different than my guess.

[Michele Selvaggi]

What you obtain is expected since V.Mag() is by construction equal to sqrt(px*px+py*py+pz*pz).
You won't be able to reconstruct the Z mass, since the missing Energy gives you only access to the SUM vector of the two neutrinos.

[Li Huang]

Hi,

sorry I can't get your point. For example, we can construct the z invariant mass by leptons, like z > l+ l- , if the final state momentum is p1 p2, then the invariant mass is p^2=(p1+p2)2 . So, if the missing energy give me the sum vector, I should get the invariant mass. I guess you mean missing Energy branch contains only px py pz ( just a 3 dimenonsion vector) , no E information. Then you calculate E = sqrt(px*px+py*py+pz*pz), not by E = init-E - detected-E ? If you do so, it is strange since the total energy may not equal to the init energy.

So, if I want to reconstruct the Z mass , what can I do , Can you help me?

[Michele Selvaggi]

Hi,

you are right the term missing energy is misleading, it should rather be called missing momentum. In your case I would simply compute as

E = sqrt(s) - Evisible

Evisible = jet[0]->E + jet[1]->E (I assume Higgs goes always to bb)

Then build a lorentz vector Vz.SetPxPyPzE(me.Px(),me.Py(),me.Pz(),E)
and Mz = Vz.M().

[Li Huang]

Hi,

Thanks for your answer! For the ideal case , I agree with Evisible = jet[0]->E + jet[1]->E (I assume Higgs goes always to bb). Since we may get a more complex result after parton shower , so I guess Evisible should sum all the visible constituents' ( in final state) energy, then, what should I do?
In the delphes result, I see that we have particle, GenJet, Track,Tower,ChargedHadron,NeutrinoHadron,Photon...etc. Shall I go through all the branches to get the visible energy?

[Michele Selvaggi]

Hi,

If you are doing the analysis with reconstructed particles (which I assume you are since you are using MissingET) you just need EFlowPhoton, EFlowNeutralHadron, and EFlowTracks. Nothing else, otherwise you double count particles.

Conversely If you want to do the analysis at generated level you can take all status = 1 gen particles.

However I just realised that some of these collections were missing in the ILD card.
I have just pushed the modification to the card, you can get the new version of Delphes with this command:

git clone https://github.com/delphes/delphes.git

Hope this helps,
Michele

[Li Huang]

Hi,

Yes,I am reconstructing particles.
Honestly, I didn't know EFlow... , I think the primary objects are : Track, Tower, EFlowTrack, EFlowPhoton, EFlowNeutralHardon, I don't know why you say nothing else.
Or I can just consider the final objects like: MissingET, ScalarHT, Jet, Muon, Electron, Photon, using these information to calculate Missing momentum myself?
Or I can just consider the generate level by Particle and GenJet ?

Maybe this question is too naive but I didn't find a manual to explain EFlow-objects and how to operate them, so it will be very nice for you to share a manual or introduction,(I thought the manual of Delphes or FastJet will explain these object: EFlowPhoton, EFlowTrack..., but it doesn't).

And I tried the new version and it contains the new branches.

[Michele Selvaggi]

EFlowObjects contain all charged particles (e,mu,charged hadrons) and neutral deposits (neutral hadrons, photons) all you need.
The missing ET is already calculated using EFLowObjects, you can actually use that as cross check.

you can read about them in our paper:

​http://arxiv.org/abs/1307.6346

The reason why you should not use directly photons, Muons etc, is that those collections have some cuts already applied.

[Li Huang]

Hi,

I am so sorry that when I want to check your last comment, there is a bug,the code is :

/*
This macro shows how to access the particle-level reference for reconstructed objects.
It is also shown how to loop over the jet constituents.

root -l examples/Example3.C'("delphes_output.root")'
*/

#ifdef CLING
RLOAD_LIBRARY(libDelphes)
#include "/home/ui/tools/package/delphes/classes/DelphesClasses.h"
#include "/home/ui/tools/package/delphes/external/ExRootAnalysis/ExRootTreeReader.h"
#include "/home/ui/tools/package/delphes/external/ExRootAnalysis/ExRootResult.h"
#else
class ExRootTreeReader;
class ExRootResult;
#endif

#include "/home/ui/tools/package/delphes/modules/SimpleCalorimeter.h"

------------------------------------------------------------------------------
class EFlowTrack;
class EFlowPhoton;
class EFlowNeutralHadron;

#include<vector>

void AnalyseEvents(ExRootTreeReader *treeReader)
{

TClonesArray *branchParticle = treeReader->UseBranch("Particle");
TClonesArray *branchElectron = treeReader->UseBranch("Electron");
TClonesArray *branchPhoton = treeReader->UseBranch("Photon");
TClonesArray *branchMuon = treeReader->UseBranch("Muon");

TClonesArray *branchTrack = treeReader->UseBranch("Track");
TClonesArray *branchTower = treeReader->UseBranch("Tower");

TClonesArray *branchEFlowTrack = treeReader->UseBranch("EFlowTrack");
TClonesArray *branchEFlowPhoton = treeReader->UseBranch("EFlowPhoton");
TClonesArray *branchEFlowNeutralHadron = treeReader->UseBranch("EFlowNeutralHadron");
TClonesArray *branchJet = treeReader->UseBranch("Jet");

TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");

Long64_t allEntries = treeReader->GetEntries();

cout << " Chain contains " << allEntries << " events" << endl;

GenParticle *particle;
Electron *electron;
Photon *photon;
Muon *muon;

Track *track;
Tower *tower;

Jet *jet;
TObject *object;

TLorentzVector momentum;

Float_t Eem, Ehad;
Bool_t skip;

Long64_t entry;

Int_t i, j, pdgCode;

EFlowTrack *eflowtrack;
EFlowPhoton *eflowphoton;
EFlowNeutralHadron *eflowneutralhadron;
TLorentzVector initmomentum,missingmomentum,finalmomentum,minusmomentum;
vector <TLorentzVector> vecmomentum;
MissingET *met;

initmomentum.SetPxPyPzE(0.0,0.0,0.0,500);
Loop over all events
for(entry = 0; entry < allEntries; ++entry)
{

Load selected branches with data from specified event
treeReader->ReadEntry(entry);
for(i = 0; i < branchEFlowTrack->GetEntriesFast(); ++i)
{

eflowtrack = (EFlowTrack*) branchEFlowTrack->At(i);
momentum.SetPtEtaPhiM(eflowtrack->PT,eflowtrack->Eta,eflowtrack->Phi,0.0);
vecmomentum.push_back(momentum);

}
for(i = 0; i < branchEFlowPhoton->GetEntriesFast(); ++i)
{

eflowphoton = (EFlowPhoton*) branchEFlowPhoton->At(i);
momentum.SetPtEtaPhiM(eflowphoton->PT,eflowphoton->Eta,eflowphoton->Phi,0.0);
vecmomentum.push_back(momentum);

}
for(i = 0; i < branchEFlowNeutralHadron->GetEntriesFast(); ++i)
{

eflowneutralhadron = (EFlowNeutralHadron*) branchEFlowNeutralHadron->At(i);
momentum.SetPtEtaPhiM(eflowneutralhadron->PT,eflowneutralhadron->Eta,eflowneutralhadron->Phi,0.0);
vecmomentum.push_back(momentum);

}
for(i = 0; i < vecmomentum.size(); i++)
{

finalmomentum = initmomentum - vecmomentum[i];

}

if(branchMissingET->GetEntriesFast() == 1 )
{

met = (MissingET*) branchMissingET->At(0);
missingmomentum.SetPtEtaPhiM(met->MET,met->Eta,met->Phi,0.0);

}

minusmomentum = finalmomentum - missingmomentum;

cout <<entry << " "<< minusmomentum.Pt() << " " << minusmomentum.Eta() << " " << minusmomentum.Phi() << endl;

}

}

------------------------------------------------------------------------------

void testecal()
{

gSystem->Load("/home/ui/tools/package/delphes/libDelphes");

TChain *chain = new TChain("Delphes");
chain->Add("/home/ui/tools/package/delphes/zh.root");

ExRootTreeReader *treeReader = new ExRootTreeReader(chain);

AnalyseEvents(treeReader);

cout << " Exiting..." << endl;

}

------------------------------------------------------------------------------

There is a error , said,
ui@ui-Lenovo-Product:~/work/eezhh/dataanalysis/finalstate$ root -l testecal.C
root [0]
Processing testecal.C...
Chain contains 50000 events
Error: Symbol eflowtrack is not defined in current scope testecal.C:84:
Error: Failed to evaluate eflowtrack->PT
Error: Symbol eflowtrack is not defined in current scope testecal.C:84:
Error: Failed to evaluate eflowtrack->Eta
Error: Symbol eflowtrack is not defined in current scope testecal.C:84:
Error: Failed to evaluate eflowtrack->Phi
* Interpreter error recovered *
root [1]

I find that I can not use the methods of the EFlowTrack,EFlowPhoton and EFlowNeutralHadron, I can declare and assign "eflowtrack" but can not use "eflowtrack->PT" ...

[Michele Selvaggi]

Hi,

sorry for the late reply.
eflowTrack is an object of type Track, not EFlowTrack. Similarly EFlowNeutralHadron and EFlowPhoton are of type Tower.

Try this instead:

Track *eflowtrack;
Tower *eflowphoton;
Tower *eflowneutralhadron;

…

eflowtrack = (Track*) branchEFlowTrack->At(i);
eflowphoton = (Tower*) branchEFlowPhoton->At(i);
eflowneutralhadron = (Tower*) branchEFlowNeutralHadron->At(i);

Michele

[Li Huang]

Hi,

Thanks a lot anyway. I tried it for another way, I just modify the source code and add a member missing mass of MissingET . I will try to use your method later.

Best,
Li