All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
Utilities.cc
Go to the documentation of this file.
1 // *******************************************************
2 // some useful functions
3 // *******************************************************
4 #include "Utilities.h"
5 #include "TVector3.h"
6 #include "TMath.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TLorentzVector.h"
10 #include <UTIL/LCRelationNavigator.h>
11 #include <IMPL/ReconstructedParticleImpl.h>
12 #include <iomanip>
13 
14 namespace isolep{
15 
16 Int_t getMCSerial(MCParticle *mcPart, LCCollection *colMCP) {
17  // get the serial code of the particle in the MC particle list
18  Int_t nMCP = colMCP->getNumberOfElements();
19  Int_t serialID;
20  for (Int_t i=0;i<nMCP;i++) {
21  MCParticle *mcp = dynamic_cast<MCParticle*>(colMCP->getElementAt(i));
22  if (mcp == mcPart) {
23  serialID = i;
24  return serialID;
25  }
26  }
27  return -1;
28 }
29 
30  MCParticle *getMCParticle(ReconstructedParticle *recPart, LCCollection *colMCTL) {
31  Double_t weight;
32  Int_t nMCTL;
33  return getMCParticle(recPart,colMCTL,weight,nMCTL);
34  }
35 
36  MCParticle *getMCParticle(ReconstructedParticle *recPart, LCCollection *colMCTL, Double_t &weight) {
37  Int_t nMCTL;
38  return getMCParticle(recPart,colMCTL,weight,nMCTL);
39  }
40 
41  MCParticle *getMCParticle(ReconstructedParticle *recPart, LCCollection *colMCTL, Double_t &weight, Int_t &nMCTL) {
42  // get the corresponding MC particle of one reconstructed particle using the MCTruthLinker information
43  MCParticle *mcLinkedParticle = NULL;
44  // Int_t iLink = -1;
45  LCRelationNavigator *navMCTL = new LCRelationNavigator(colMCTL);
46  LCObjectVec vecMCTL = navMCTL->getRelatedToObjects(recPart);
47  FloatVec vecWgtMCTL = navMCTL->getRelatedToWeights(recPart);
48  weight = 0.;
49  nMCTL = vecMCTL.size();
50  for (Int_t i=0;i<nMCTL;i++) {
51  MCParticle *mcPart = dynamic_cast<MCParticle *>(vecMCTL[i]);
52  if (vecWgtMCTL[i] > weight) { // find the linked particle with largest weight as the mc truth particle
53  weight = vecWgtMCTL[i];
54  mcLinkedParticle = mcPart;
55  }
56  }
57  delete navMCTL;
58  return mcLinkedParticle;
59  }
60 
61 Int_t getLinkedMCParticle(ReconstructedParticle *recPart, LCCollection *colMCTL, Double_t &weight, Int_t &nMCTL) {
62  // get the corresponding MC particle of one reconstructed particle using the MCTruthLinker information
63  // MCParticle *mcLinkedParticle = NULL;
64  Int_t iLink = -1;
65  LCRelationNavigator *navMCTL = new LCRelationNavigator(colMCTL);
66  LCObjectVec vecMCTL = navMCTL->getRelatedToObjects(recPart);
67  FloatVec vecWgtMCTL = navMCTL->getRelatedToWeights(recPart);
68  Double_t mcEnergyMax = -1.0;
69  weight = 0.;
70  nMCTL = vecMCTL.size();
71  for (Int_t i=0;i<nMCTL;i++) {
72  MCParticle *mcPart = dynamic_cast<MCParticle *>(vecMCTL[i]);
73  Double_t mcEnergy = mcPart->getEnergy();
74  if (mcEnergy > mcEnergyMax) { // find the linked particle with largest energy as the mc truth particle
75  mcEnergyMax = mcEnergy;
76  // mcLinkedParticle = mcPart;
77  weight = vecWgtMCTL[i];
78  iLink = i;
79  }
80  }
81  // return mcLinkedParticle;
82  return iLink;
83 }
84 
85 Int_t getOriginalPDG(MCParticle *mcPart, Bool_t iHiggs) {
86  // get the PDG of the original particle where the MCParticle comes from
87  // Int_t nParents = mcPart->getNumberOfParents();
88  Int_t nParents = mcPart->getParents().size();
89  while (nParents > 0) {
90  // MCParticle *mother = mcPart->getParent(0);
91  MCParticle *mother = mcPart->getParents()[0];
92  // nParents = mother->getNumberOfParents();
93  nParents = mother->getParents().size();
94  mcPart = mother;
95  Int_t pdg = mcPart->getPDG();
96  if (!iHiggs) {
97  if (abs(pdg) == 24 || abs(pdg) == 23 || abs(pdg) == 25) break;
98  }
99  else {
100  if (abs(pdg) == 25) break;
101  }
102  }
103  Int_t originalPDG = mcPart->getPDG();
104  return originalPDG;
105 }
106 
107  Int_t getOriginalPDGForIsoLep(MCParticle *mcPart) {
108  // get the PDG of the original particle where the MCParticle comes from
109  // Int_t nParents = mcPart->getNumberOfParents();
110  Int_t nParents = mcPart->getParents().size();
111  // if (mcPart->isOverlay()) return -22;
112  Double_t charge = mcPart->getCharge();
113  while (nParents > 0) {
114  MCParticle *mother = mcPart->getParents()[0];
115  if (nParents > 1) {
116  for (Int_t i=0;i<nParents;i++) {
117  MCParticle *temp = mcPart->getParents()[i];
118  if (TMath::Abs(temp->getCharge() - charge) < 0.1) {
119  mother = temp;
120  break;
121  }
122  }
123  }
124  mcPart = mother;
125  nParents = mother->getParents().size();
126  Int_t pdg = mcPart->getPDG();
127  if (abs(pdg) == 13 || abs(pdg) == 11 || abs(pdg) == 15) {
128  if (mcPart->getParents().size() > 0) {
129  MCParticle *mmother = mcPart->getParents()[0];
130  Int_t mpdg = mmother->getPDG();
131  if (abs(mpdg) == 11) {
132  // std::cerr << "Mother is an electron!" << endl;
133  break;
134  }
135  }
136  }
137  else if (abs(pdg) == 24 || abs(pdg) == 23 || pdg == 25 || pdg == 21 || abs(pdg) == 15) {
138  break;
139  }
140  else if (abs(pdg) <= 6 && abs(pdg) >= 1) {
141  break;
142  }
143 
144  }
145  // if (mcPart->isOverlay()) return -22;
146  Int_t originalPDG = mcPart->getPDG();
147  return originalPDG;
148  }
149 
150  Int_t getOriginalPDGForIsoLep(MCParticle *mcPart, LCCollection *colMC) {
151  // get the PDG of the original particle where the MCParticle comes from
152  Int_t nParents = mcPart->getParents().size();
153  Int_t serial = getMCSerial(mcPart,colMC);
154  if (mcPart->isOverlay()) return -22;
155  MCParticle *previous = mcPart;
156  while (nParents > 0 && serial > 11) {
157  Double_t charge = previous->getCharge();
158  MCParticle *mother = mcPart->getParents()[0];
159  if (nParents > 1) {
160  for (Int_t i=0;i<nParents;i++) {
161  MCParticle *temp = mcPart->getParents()[i];
162  if (TMath::Abs(temp->getCharge() - charge) < 0.1) {
163  mother = temp;
164  break;
165  }
166  }
167  }
168  previous = mcPart;
169  mcPart = mother;
170  if (mcPart->isOverlay()) break;
171  nParents = mcPart->getParents().size();
172  serial = getMCSerial(mcPart,colMC);
173  }
174  if (mcPart->isOverlay()) return -22;
175  Int_t originalPDG = mcPart->getPDG();
176  return originalPDG;
177  }
178 
179 Int_t getOriginalSerial(MCParticle *mcPart, LCCollection *colMCP, Bool_t iHiggs) {
180  // get the serial number of the original particle where the MCParticle comes from
181  // Int_t nParents = mcPart->getNumberOfParents();
182  Int_t nParents = mcPart->getParents().size();
183  while (nParents > 0) {
184  // MCParticle *mother = mcPart->getParent(0);
185  MCParticle *mother = mcPart->getParents()[0];
186  // nParents = mother->getNumberOfParents();
187  nParents = mother->getParents().size();
188  mcPart = mother;
189  Int_t pdg = mcPart->getPDG();
190  if (!iHiggs) {
191  if (abs(pdg) == 24 || abs(pdg) == 23 || abs(pdg) == 25) break;
192  }
193  else {
194  if (abs(pdg) == 25) break;
195  }
196  }
197  Int_t originalSerial = getMCSerial(mcPart,colMCP);
198  return originalSerial;
199 }
200 
201 Int_t getOriginalSerial(ReconstructedParticle *recPart, LCCollection *colMCTL, LCCollection *colMCP, Bool_t iHiggs) {
202  // get the serial number of the original particle where the PFO comes from
203 
204  Int_t originalSerial = -1;
205  LCRelationNavigator *navMCTL = new LCRelationNavigator(colMCTL);
206  LCObjectVec vecMCTL = navMCTL->getRelatedToObjects(recPart);
207  if (vecMCTL.size() > 0) {
208  MCParticle *mcPart = dynamic_cast<MCParticle *>(vecMCTL[0]);
209  originalSerial = getOriginalSerialForZHH(mcPart,colMCP);
210  }
211  return originalSerial;
212 }
213 
214 Int_t getOriginalSerialForZHH(MCParticle *mcPart, LCCollection *colMCP) {
215  // get the serial number of the original particle where the MCParticle comes from
216  // Int_t nParents = mcPart->getNumberOfParents();
217  Int_t nParents = mcPart->getParents().size();
218  while (nParents > 0) {
219  // MCParticle *mother = mcPart->getParent(0);
220  MCParticle *mother = mcPart->getParents()[0];
221  // nParents = mother->getNumberOfParents();
222  nParents = mother->getParents().size();
223  mcPart = mother;
224  Int_t pdg = mcPart->getPDG();
225  if (abs(pdg) == 25) {
226  if (mcPart->getParents().size() > 0) {
227  MCParticle *mmother = mcPart->getParents()[0];
228  Int_t mpdg = mmother->getPDG();
229  if (mpdg == 11) break;
230  }
231  }
232  }
233  Int_t originalSerial = getMCSerial(mcPart,colMCP);
234  return originalSerial;
235 }
236 
237 Int_t getLeptonID(ReconstructedParticle *recPart) {
238  // electron identification using ratios of energies deposited in ECal, HCal and Momentum
239  Int_t iLeptonType = 0;
240  Double_t fElectronCut1 = 0.5; // lower edge of totalCalEnergy/momentum
241  Double_t fElectronCut2 = 1.3; // upper edge of totalCalEnergy/momentum
242  Double_t fElectronCut3 = 0.9; // lower edge of ecalEnergy/totalCalEnergy
243  Double_t fMuonCut1 = 0.3; // upper edge of totalCalEnergy/momentum
244  // Double_t fMuonCut2 = 0.5; // upper edge of ecalEnergy/totalCalEnergy
245  Double_t fMuonCut3 = 1.2; // lower edge of yoke energy
246  Double_t fPhotonCut1 = 0.7; // lower edge of totalCalEnergy/momentum
247  Double_t fPhotonCut2 = 1.3; // upper edge of totalCalEnergy/momentum
248  Double_t fPhotonCut3 = 0.9; // lower edge of ecalEnergy/totalCalEnergy
249 
250  // Double_t energy = recPart->getEnergy();
251  Double_t charge = recPart->getCharge();
252  Double_t ecalEnergy = 0;
253  Double_t hcalEnergy = 0;
254  Double_t yokeEnergy = 0;
255  Double_t totalCalEnergy = 0;
256  // Int_t nHits = 0;
257  std::vector<lcio::Cluster*> clusters = recPart->getClusters();
258  for (std::vector<lcio::Cluster*>::const_iterator iCluster=clusters.begin();iCluster!=clusters.end();++iCluster) {
259  ecalEnergy += (*iCluster)->getSubdetectorEnergies()[0];
260  hcalEnergy += (*iCluster)->getSubdetectorEnergies()[1];
261  yokeEnergy += (*iCluster)->getSubdetectorEnergies()[2];
262  ecalEnergy += (*iCluster)->getSubdetectorEnergies()[3];
263  hcalEnergy += (*iCluster)->getSubdetectorEnergies()[4];
264  // CalorimeterHitVec calHits = (*iCluster)->getCalorimeterHits();
265  // nHits = calHits.size();
266  }
267  totalCalEnergy = ecalEnergy + hcalEnergy;
268  TVector3 momentum = TVector3(recPart->getMomentum());
269  Double_t momentumMagnitude = momentum.Mag();
270  Double_t fEpsilon = 1.E-10;
271  if (totalCalEnergy/momentumMagnitude > fElectronCut1 && totalCalEnergy/momentumMagnitude < fElectronCut2 &&
272  ecalEnergy/(totalCalEnergy + fEpsilon) > fElectronCut3 && TMath::Abs(charge) > 0.5) {
273  iLeptonType = 11;
274  }
275  else if (totalCalEnergy/momentumMagnitude > fPhotonCut1 && totalCalEnergy/momentumMagnitude < fPhotonCut2 &&
276  ecalEnergy/(totalCalEnergy + fEpsilon) > fPhotonCut3 && TMath::Abs(charge) < 0.5) {
277  iLeptonType = 22;
278  }
279  // else if (TMath::Abs(charge) > 0.5 && totalCalEnergy/momentumMagnitude < fMuonCut1 &&
280  // ecalEnergy/(totalCalEnergy + fEpsilon) < fMuonCut2) {
281  else if (TMath::Abs(charge) > 0.5 && totalCalEnergy/momentumMagnitude < fMuonCut1 &&
282  yokeEnergy > fMuonCut3) {
283  iLeptonType = 13;
284  }
285 
286  return iLeptonType;
287 }
288 
289 Bool_t getFSRTag(ReconstructedParticle *motherPart, ReconstructedParticle *recPart, Double_t fCosFSRCut) {
290  // tag the particle recPart if it is from the Bremmstrahlung or Final-State-Radiation of another particle motherPart
291  Bool_t isFSR = kFALSE;
292  // Double_t charge = motherPart->getCharge(); // mother particle should be charged
293  TVector3 momentumMother = TVector3(motherPart->getMomentum());
294  TVector3 momentum = TVector3(recPart->getMomentum());
295  Double_t cosFSR = momentum.Dot(momentumMother)/momentum.Mag()/momentumMother.Mag();
296  Int_t iType = getLeptonID(recPart);
297  Int_t iTypeMother = getLeptonID(motherPart);
298  // if (TMath::Abs(charge) > 0.5 && cosFSR > fCosFSRCut) {
299  // if (iType == 11 || iType == 22) {
300  if (cosFSR > fCosFSRCut && (iTypeMother == 11 || iTypeMother == 13)) {
301  if (iType == 22) {
302  isFSR = kTRUE;
303  }
304  }
305  return isFSR;
306 }
307 
308 Bool_t getSplitTag(ReconstructedParticle *motherPart, ReconstructedParticle *recPart) {
309  // developed from Mark Tomthon's ZFinder Processor
310  // tag the particle recPart if it is the cluster splited from motherPart
311  Bool_t isSplit = kFALSE;
312 
313  // information of mother particle
314  Double_t trackMom=0.;
315  Double_t clusterMom=0.;
316  Double_t sigmaMom=999.;
317  Double_t chiMom=999.;
318  Double_t sigpMom=0.;
319  TVector3 vecMom;
320  const EVENT::ClusterVec ci = motherPart->getClusters();
321  const EVENT::TrackVec ti = motherPart->getTracks();
322  trackMom = motherPart->getEnergy();
323  if(ti.size()==1)sigpMom = trackMom*sqrt(ti[0]->getCovMatrix()[5])/TMath::Abs(ti[0]->getOmega());
324  if(ci.size()>0){
325  clusterMom = ci[0]->getEnergy();
326  sigmaMom = 0.18*sqrt(clusterMom);
327  chiMom = (trackMom-clusterMom)/TMath::Sqrt(sigmaMom*sigmaMom+sigpMom*sigpMom);
328  vecMom = TVector3(ci[0]->getPosition()[0],ci[0]->getPosition()[1],ci[0]->getPosition()[2]);
329  }
330 
331  // calculate the distance between cluster and the mother particle cluster
332  Double_t dr = 999.;
333  const EVENT::ClusterVec c = recPart->getClusters();
334  if(c.size()==1){
335  TVector3 vecg(c[0]->getPosition()[0],c[0]->getPosition()[1],c[0]->getPosition()[2]);
336  TVector3 v = vecg.Cross(vecMom);
337  float magg = vecg.Mag();
338  if(magg>0) dr = v.Mag()/magg;
339  }
340 
341  // criteria for tag
342  Double_t eg = recPart->getEnergy();
343  float chiNew = (trackMom-clusterMom-eg)/sigmaMom;
344  if(dr<20.0) isSplit = true;
345  // if fairly close merge if chi2 for matching improves greatly
346  if(dr<30.0 && chiMom>4.0 && TMath::Abs(chiNew)<chiMom) isSplit = true;
347  if(dr<40.0 && chiMom>5.0 && TMath::Abs(chiNew)<chiMom) isSplit = true;
348  if(dr<50.0 && chiMom>7.0 && TMath::Abs(chiNew)<chiMom) isSplit = true;
349  // sanity check
350  if(TMath::Abs(chiMom)<2.0 && chiNew*chiNew>chiMom*chiMom+5.0) isSplit = false;
351  // always merge if very close - can't expect reconstruction to work
352  if(dr<10.0) isSplit = true;
353 
354  // empirical correction
355  Double_t costheta = getCosTheta(motherPart,recPart);
356  if (costheta<0.99995 && eg/trackMom < 0.03) isSplit = true;
357 
358  return isSplit;
359 }
360 
361 Double_t getConeEnergy(ReconstructedParticle *recPart, LCCollection *colPFO, Double_t cosCone) {
362  // get the cone energy of the particle
363  Int_t nPFOs = colPFO->getNumberOfElements();
364  Double_t coneEnergy = 0.;
365  TVector3 momentum0 = TVector3(recPart->getMomentum());
366  for (Int_t i=0;i<nPFOs;i++) {
367  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
368  if (pfo != recPart) {
369  TVector3 momentum = TVector3(pfo->getMomentum());
370  Double_t cosTheta = momentum0.Dot(momentum)/momentum0.Mag()/momentum.Mag();
371  if (cosTheta > cosCone) {
372  coneEnergy += pfo->getEnergy();
373  }
374  }
375  }
376  return coneEnergy;
377 }
378 
379 void getConeEnergy(ReconstructedParticle *recPart, LCCollection *colPFO, Double_t cosCone, Bool_t woFSR, Double_t coneEnergy[3], Double_t pFSR[4]) {
380  // get the cone energy of the particle
381  // woFSR = kTRUE;
382  TLorentzVector lortzFSR = TLorentzVector(0, 0, 0, 0);
383  Int_t nPFOs = colPFO->getNumberOfElements();
384  TVector3 momentum0 = TVector3(recPart->getMomentum());
385  for (Int_t i=0;i<nPFOs;i++) {
386  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
387  if (pfo == recPart) continue;
388  Double_t energy = pfo->getEnergy();
389  TVector3 momentum = TVector3(pfo->getMomentum());
390  if (woFSR) {
391  Bool_t isFSR = getFSRTag(recPart,pfo);
392  if (isFSR) {
393  lortzFSR += TLorentzVector(momentum,energy);
394  continue;
395  }
396  }
397  Double_t charge = pfo->getCharge();
398  Double_t cosTheta = momentum0.Dot(momentum)/momentum0.Mag()/momentum.Mag();
399  if (cosTheta > cosCone) {
400  coneEnergy[0] += energy;
401  if (TMath::Abs(charge) < 0.5) {
402  coneEnergy[1] += energy;
403  }
404  else {
405  coneEnergy[2] += energy;
406  }
407  }
408  }
409  pFSR[0] = lortzFSR.Px();
410  pFSR[1] = lortzFSR.Py();
411  pFSR[2] = lortzFSR.Pz();
412  pFSR[3] = lortzFSR.E();
413 
414  return;
415 }
416 
417 void getConeEnergy(ReconstructedParticle *recPart, LCCollection *colPFO, Double_t cosCone, Bool_t woFSR, Double_t coneEnergy[3], Double_t pFSR[4],
418  Double_t cosCone2, Double_t pCone2[4], Int_t &nConePhoton ) {
419  // get the cone energy of the particle
420  // add another larger cone
421  // woFSR = kTRUE;
422  TLorentzVector lortzFSR = TLorentzVector(0, 0, 0, 0);
423  TLorentzVector lortzCon = TLorentzVector(0, 0, 0, 0);
424  Int_t nPFOs = colPFO->getNumberOfElements();
425  TVector3 momentum0 = TVector3(recPart->getMomentum());
426  nConePhoton = 0;
427  for (Int_t i=0;i<nPFOs;i++) {
428  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
429  if (pfo == recPart) continue;
430  Double_t energy = pfo->getEnergy();
431  TVector3 momentum = TVector3(pfo->getMomentum());
432  if (woFSR) {
433  Bool_t isFSR = getFSRTag(recPart,pfo);
434  if (isFSR) {
435  lortzFSR += TLorentzVector(momentum,energy);
436  Int_t iType = getLeptonID(pfo);
437  if (iType == 22) nConePhoton++;
438  continue;
439  }
440  }
441  Double_t charge = pfo->getCharge();
442  Double_t cosTheta = momentum0.Dot(momentum)/momentum0.Mag()/momentum.Mag();
443  if (cosTheta > cosCone) {
444  coneEnergy[0] += energy;
445  if (TMath::Abs(charge) < 0.5) {
446  coneEnergy[1] += energy;
447  }
448  else {
449  coneEnergy[2] += energy;
450  }
451  Int_t iType = getLeptonID(pfo);
452  if (iType == 22) nConePhoton++;
453  }
454  if (cosTheta > cosCone2) {
455  lortzCon += TLorentzVector(momentum,energy);
456  }
457  }
458  pFSR[0] = lortzFSR.Px();
459  pFSR[1] = lortzFSR.Py();
460  pFSR[2] = lortzFSR.Pz();
461  pFSR[3] = lortzFSR.E();
462  pCone2[0] = lortzCon.Px();
463  pCone2[1] = lortzCon.Py();
464  pCone2[2] = lortzCon.Pz();
465  pCone2[3] = lortzCon.E();
466  return;
467 }
468 
469 TLorentzVector getFSRMomentum(ReconstructedParticle *recPart, LCCollection *colPFO) {
470  // get the cone energy of the particle
471  TLorentzVector lortzFSR = TLorentzVector(0.,0.,0.,0.);
472  Int_t nPFOs = colPFO->getNumberOfElements();
473  TVector3 momentum0 = TVector3(recPart->getMomentum());
474  for (Int_t i=0;i<nPFOs;i++) {
475  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
476  if (pfo == recPart) continue;
477  Bool_t isFSR = getFSRTag(recPart,pfo);
478  if (! isFSR) continue;
479  Double_t energy = pfo->getEnergy();
480  TVector3 momentum = TVector3(pfo->getMomentum());
481  lortzFSR += TLorentzVector(momentum,energy);
482  }
483 
484  return lortzFSR;
485 }
486 
487 Double_t getConeEnergy(ReconstructedParticle *recPart, LCCollection *colPFO, Double_t cosCone, Int_t mode) {
488  // get the cone energy of the particle
489  Int_t nPFOs = colPFO->getNumberOfElements();
490  Double_t coneEnergy = 0.,coneEnergyC = 0.,coneEnergyN = 0.;
491  TVector3 momentum0 = TVector3(recPart->getMomentum());
492  for (Int_t i=0;i<nPFOs;i++) {
493  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
494  if (pfo != recPart) {
495  TVector3 momentum = TVector3(pfo->getMomentum());
496  Int_t iCharge = pfo->getCharge();
497  Double_t cosTheta = momentum0.Dot(momentum)/momentum0.Mag()/momentum.Mag();
498  if (cosTheta > cosCone) {
499  coneEnergy += pfo->getEnergy();
500  if (iCharge == 0) {
501  coneEnergyN += pfo->getEnergy();
502  }
503  else {
504  coneEnergyC += pfo->getEnergy();
505  }
506  }
507  }
508  }
509  if (mode == 0) {
510  return coneEnergy;
511  }
512  else if (mode == 1) {
513  return coneEnergyC;
514  }
515  else if (mode == 2) {
516  return coneEnergyN;
517  }
518  else {
519  return 99999.;
520  }
521 }
522 
523 Double_t getConeEnergy(ReconstructedParticle *recPart, LCCollection *colPFO, Double_t cosCone,
524  std::vector<lcio::ReconstructedParticle*> &conePFOs) {
525  // get the cone energy of the particle
526  Int_t nPFOs = colPFO->getNumberOfElements();
527  Double_t coneEnergy = 0.;
528  TVector3 momentum0 = TVector3(recPart->getMomentum());
529  for (Int_t i=0;i<nPFOs;i++) {
530  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
531  if (pfo != recPart) {
532  TVector3 momentum = TVector3(pfo->getMomentum());
533  Double_t cosTheta = momentum0.Dot(momentum)/momentum0.Mag()/momentum.Mag();
534  if (cosTheta > cosCone) {
535  coneEnergy += pfo->getEnergy();
536  conePFOs.push_back(pfo);
537  }
538  }
539  }
540  return coneEnergy;
541 }
542 
543 void getConeEnergy(ReconstructedParticle *recPart, LCCollection *colPFO, Double_t cosCone,
544  Double_t coneEnergy[3], Double_t cosCone2, Double_t pCone2[4]){
545  // get the cone energy of the particle
546  // add another larger cone
547  TLorentzVector lortzCon = TLorentzVector(0, 0, 0, 0);
548  Int_t nPFOs = colPFO->getNumberOfElements();
549  TVector3 momentum0 = TVector3(recPart->getMomentum());
550  for (Int_t i=0;i<nPFOs;i++) {
551  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
552  if (pfo == recPart) continue;
553  Double_t energy = pfo->getEnergy();
554  TVector3 momentum = TVector3(pfo->getMomentum());
555  Double_t charge = pfo->getCharge();
556  Double_t cosTheta = momentum0.Dot(momentum)/momentum0.Mag()/momentum.Mag();
557  if (cosTheta > cosCone) {
558  coneEnergy[0] += energy;
559  if (TMath::Abs(charge) < 0.5) {
560  coneEnergy[1] += energy;
561  }
562  else {
563  coneEnergy[2] += energy;
564  }
565  }
566  if (cosTheta > cosCone2) {
567  lortzCon += TLorentzVector(momentum,energy);
568  }
569  }
570  pCone2[0] = lortzCon.Px();
571  pCone2[1] = lortzCon.Py();
572  pCone2[2] = lortzCon.Pz();
573  pCone2[3] = lortzCon.E();
574  return;
575 }
576 
577 Double_t getInvariantMass(ReconstructedParticle *recPart1, ReconstructedParticle *recPart2) {
578  // get the invariant mass of two particles
579  TVector3 momentum1 = TVector3(recPart1->getMomentum());
580  TVector3 momentum2 = TVector3(recPart2->getMomentum());
581  Double_t energy1 = recPart1->getEnergy();
582  Double_t energy2 = recPart2->getEnergy();
583  Double_t invariantMass = 0.;
584  Double_t invariantMass2 = (energy1+energy2)*(energy1+energy2)-(momentum1+momentum2).Mag2();
585  if (invariantMass2 > 0.) {
586  invariantMass = sqrt(invariantMass2);
587  }
588  else {
589  invariantMass = -sqrt(TMath::Abs(invariantMass2));
590  }
591  return invariantMass;
592 }
593 
594 Double_t getInvariantMass(ReconstructedParticle *recPart1, ReconstructedParticle *recPart2, ReconstructedParticle *recPart3) {
595  // get the invariant mass of three particles
596  TVector3 momentum1 = TVector3(recPart1->getMomentum());
597  TVector3 momentum2 = TVector3(recPart2->getMomentum());
598  TVector3 momentum3 = TVector3(recPart3->getMomentum());
599  Double_t energy1 = recPart1->getEnergy();
600  Double_t energy2 = recPart2->getEnergy();
601  Double_t energy3 = recPart3->getEnergy();
602  Double_t invariantMass = 0.;
603  Double_t invariantMass2 = (energy1+energy2+energy3)*(energy1+energy2+energy3)-(momentum1+momentum2+momentum3).Mag2();
604  if (invariantMass2 > 0.) {
605  invariantMass = sqrt(invariantMass2);
606  }
607  else {
608  invariantMass = -sqrt(TMath::Abs(invariantMass2));
609  }
610  return invariantMass;
611 }
612 
613 Int_t isSelectedByFastJet( ReconstructedParticle *pfo, LCCollection *colFastJet, Double_t &ratioEPartEJet, Double_t &ratioPTMJet ) {
614  // check the PFO if it is selectred by FastJet clustering
615  Int_t iFastJet = 0;
616  Int_t nJets = colFastJet->getNumberOfElements();
617  Int_t iJet = -1;
618  for (Int_t i=0;i<nJets;i++) {
619  ReconstructedParticle *jet = dynamic_cast<ReconstructedParticle*>(colFastJet->getElementAt(i));
620  std::vector<lcio::ReconstructedParticle*> partVec = jet->getParticles();
621  for (std::vector<lcio::ReconstructedParticle*>::const_iterator iPart=partVec.begin();iPart!=partVec.end();++iPart) {
622  if ((*iPart) == pfo) {
623  iFastJet = 1;
624  iJet = i;
625  break;
626  }
627  }
628  }
629  if (iJet >= 0) {
630  ReconstructedParticle *theJet = dynamic_cast<ReconstructedParticle*>(colFastJet->getElementAt(iJet)); // the jet where the pfo belongs
631  // get the variables used by LAL Lepton Finder
632  Double_t ePart = pfo->getEnergy();
633  TVector3 pPart = pfo->getMomentum();
634  Double_t eJet = theJet->getEnergy();
635  TVector3 pJet = theJet->getMomentum();
636  // Double_t mJet = theJet->getMass();
637  Double_t mJet = eJet*eJet > pJet.Mag2() ? TMath::Sqrt(eJet*eJet-pJet.Mag2()) : TMath::Sqrt(-eJet*eJet+pJet.Mag2());
638  ratioEPartEJet = ePart/eJet;
639  ratioPTMJet = pPart.Pt(pJet)/mJet;
640  }
641 
642  return iFastJet;
643 }
644 
645 void doPhotonRecovery(ReconstructedParticle *electron, LCCollection *colPFO, ReconstructedParticleImpl *recoElectron, Double_t fCosFSRCut) {
646  // recover the BS and FSR photons
647  TLorentzVector lortzElectron = TLorentzVector(electron->getMomentum(),electron->getEnergy());
648  Int_t nPFOs = colPFO->getNumberOfElements();
649  for (Int_t i=0;i<nPFOs;i++) {
650  ReconstructedParticle *recPart = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
651  if (recPart == electron) continue;
652  Bool_t isFSR = getFSRTag(electron,recPart,fCosFSRCut);
653  if (! isFSR) continue;
654  recoElectron->addParticle(recPart);
655  Bool_t isSplit = getSplitTag(electron,recPart);
656  if (isSplit) continue;
657  lortzElectron += TLorentzVector(recPart->getMomentum(),recPart->getEnergy());
658  }
659  Double_t energy = lortzElectron.E();
660  Double_t mass = lortzElectron.M();
661  Double_t momentum[3] = {lortzElectron.Px(),lortzElectron.Py(),lortzElectron.Pz()};
662  Double_t charge = electron->getCharge();
663  recoElectron->setMomentum(momentum);
664  recoElectron->setEnergy(energy);
665  recoElectron->setMass(mass);
666  recoElectron->setCharge(charge);
667  recoElectron->setType(94);
668 
669 }
670 
671 void doPhotonRecovery(ReconstructedParticle *electron, LCCollection *colPFO, ReconstructedParticleImpl *recoElectron, Double_t fCosFSRCut,
672  Int_t lepType, std::vector<lcio::ReconstructedParticle*> &photons) {
673  // recover the BS and FSR photons
674  TLorentzVector lortzElectron = TLorentzVector(electron->getMomentum(),electron->getEnergy());
675  recoElectron->addParticle(electron);
676  Int_t nPFOs = colPFO->getNumberOfElements();
677  for (Int_t i=0;i<nPFOs;i++) {
678  ReconstructedParticle *recPart = dynamic_cast<ReconstructedParticle*>(colPFO->getElementAt(i));
679  if (recPart == electron) continue;
680  if (isFoundInVector(recPart,photons)) continue;
681  Bool_t isFSR = getFSRTag(electron,recPart,fCosFSRCut);
682  if (! isFSR) continue;
683  photons.push_back(recPart);
684  recoElectron->addParticle(recPart);
685  if (lepType == 11) {
686  // do split algorithm only for electron
687  Bool_t isSplit = getSplitTag(electron,recPart);
688  if (isSplit) continue;
689  }
690  else if (lepType == 13) {
691  }
692  lortzElectron += TLorentzVector(recPart->getMomentum(),recPart->getEnergy());
693  }
694  Double_t energy = lortzElectron.E();
695  Double_t mass = lortzElectron.M();
696  Double_t momentum[3] = {lortzElectron.Px(),lortzElectron.Py(),lortzElectron.Pz()};
697  Double_t charge = electron->getCharge();
698  recoElectron->setMomentum(momentum);
699  recoElectron->setEnergy(energy);
700  recoElectron->setMass(mass);
701  recoElectron->setCharge(charge);
702  recoElectron->setType(94);
703 
704 }
705 
706 Bool_t isFoundInVector(ReconstructedParticle *pfo, std::vector<lcio::ReconstructedParticle*> &pfos) {
707  for (std::vector<lcio::ReconstructedParticle*>::const_iterator ipfo=pfos.begin();ipfo!=pfos.end();++ipfo) {
708  if ((*ipfo) == pfo) {
709  return kTRUE;
710  }
711  }
712  return kFALSE;
713 }
714 
715 Double_t jetFunction( TLorentzVector lortz, float beta, float power) {
716  // extended Georgi Jet Function
717 
718  Double_t energy = lortz.E();
719  Double_t momentum = lortz.P();
720 
721  if (!energy) {
722  return -9999999.;
723  }
724  else {
725  return TMath::Power(energy,power)*(1.-beta+beta*momentum*momentum/energy/energy);
726  }
727 
728 }
729 
730 Double_t jetFunction( TLorentzVector lortz, float virtual_scale) {
731  // Junping Jet Function
732 
733  Double_t energy = lortz.E();
734  Double_t mass = lortz.M();
735 
736  if (!energy) {
737  return -9999999.;
738  }
739  else {
740  return energy*(1. - mass/virtual_scale);
741  }
742 
743 }
744 
745 Int_t calculateEnergyComponents(LCCollection *colMC, LCCollection *colMCTL, ReconstructedParticle *jet, Double_t energyComponents[4]) {
746  // calculate energies from different color singlets in a jet
747  // [0]: from first Higgs; [1]: from second Higgs; [2]: from other color-singlet (mostly Z); [3]: from overlay
748 
749  for (Int_t i=0;i<4;i++) {
750  energyComponents[i] = 0.;
751  }
752 
753  LCRelationNavigator *navMCTL = new LCRelationNavigator(colMCTL);
754 
755  Int_t np = jet->getParticles().size();
756  // cerr << "Debug: Npar in jet = " << np << endl;
757  for (Int_t i=0;i<np;i++) {
758  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle *>(jet->getParticles()[i]);
759  Double_t energy = pfo->getEnergy();
760  LCObjectVec vecMCTL = navMCTL->getRelatedToObjects(pfo);
761  // cerr << "Debug: energy = " << energy << " ; mclink = " << vecMCTL.size() << endl;
762  if (vecMCTL.size() > 0) {
763  MCParticle *mcPart = dynamic_cast<MCParticle *>(vecMCTL[0]);
764  Int_t originalSerial = getOriginalSerialForZHH(mcPart,colMC);
765  // cerr << "Debug: Original Serial = " << originalSerial << endl;
766  if (originalSerial == 8) { // 1st H
767  energyComponents[0] += energy;
768  }
769  else if (originalSerial == 9) { // 2nd H
770  energyComponents[1] += energy;
771  }
772  else if (mcPart->isOverlay()) { // overlay
773  energyComponents[3] += energy;
774  }
775  else {
776  energyComponents[2] += energy;
777  }
778  }
779  }
780 
781  // find the largest energy component
782  Int_t iMax = -1;
783  Double_t energyMax = -1.;
784  for (Int_t i=0;i<4;i++) {
785  if (energyComponents[i] > energyMax) {
786  energyMax = energyComponents[i];
787  iMax = i;
788  }
789  // fraction
790  // energyComponents[i] /= jet->getEnergy();
791  }
792 
793  return iMax+1;
794 }
795 
796 Int_t calculateEnergyComponents(LCCollection *colMC, LCCollection *colMCTL, ReconstructedParticle *jet, Double_t energyComponents[4],
797  Int_t iColor1, Int_t iColor2, Int_t iColor3) {
798  // calculate energies from different color singlets in a jet
799  // [0]: from first Higgs; [1]: from second Higgs; [2]: from other color-singlet (mostly Z); [3]: from overlay
800  // for qqH
801 
802  for (Int_t i=0;i<4;i++) {
803  energyComponents[i] = 0.;
804  }
805 
806  LCRelationNavigator *navMCTL = new LCRelationNavigator(colMCTL);
807 
808  Int_t np = jet->getParticles().size();
809  // cerr << "Debug: Npar in jet = " << np << endl;
810  for (Int_t i=0;i<np;i++) {
811  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle *>(jet->getParticles()[i]);
812  Double_t energy = pfo->getEnergy();
813  LCObjectVec vecMCTL = navMCTL->getRelatedToObjects(pfo);
814  // cerr << "Debug: energy = " << energy << " ; mclink = " << vecMCTL.size() << endl;
815  if (vecMCTL.size() > 0) {
816  MCParticle *mcPart = dynamic_cast<MCParticle *>(vecMCTL[0]);
817  Int_t originalSerial = getOriginalSerialForZHH(mcPart,colMC);
818  // cerr << "Debug: Original Serial = " << originalSerial << endl;
819  if (originalSerial == iColor1) { // H
820  energyComponents[0] += energy;
821  }
822  else if (originalSerial == iColor2 || originalSerial == iColor3) { // q or q-bar from Z
823  energyComponents[1] += energy;
824  }
825  else if (mcPart->isOverlay()) { // overlay
826  energyComponents[3] += energy;
827  }
828  else {
829  energyComponents[2] += energy;
830  }
831  }
832  }
833 
834  // find the largest energy component
835  Int_t iMax = -1;
836  Double_t energyMax = -1.;
837  for (Int_t i=0;i<4;i++) {
838  if (energyComponents[i] > energyMax) {
839  energyMax = energyComponents[i];
840  iMax = i;
841  }
842  // fraction
843  // energyComponents[i] /= jet->getEnergy();
844  }
845 
846  return iMax+1;
847 }
848 
849 Double_t getLikelihood(TString fname, TString hist, Double_t mass)
850 {
851 
852  TFile file(fname);
853  TH1D *hmass = dynamic_cast<TH1D *>(file.Get(hist));
854  Int_t nbin = hmass->GetNbinsX();
855  Double_t ntot = hmass->GetEntries();
856  Double_t epsilon = 0.001;
857  Double_t prob = 0.;
858  for (Int_t j=0;j<nbin;j++) {
859  Double_t nj = hmass->GetBinContent(j+1);
860  Double_t tj = hmass->GetBinCenter(j+1);
861  Double_t delta = hmass->GetBinWidth(j+1);
862  Double_t hj = TMath::Power(4./3/ntot,1./5)*delta*TMath::Sqrt(ntot/(nj+epsilon));
863  Double_t density = 1.0*nj/ntot/TMath::Sqrt(2*3.1416)/hj*TMath::Exp(-(mass-tj)*(mass-tj)/2/hj/hj);
864  // cerr << "hj: " << hj << " nj: " << nj << " tj: " << tj
865  // << " delta: " << delta << " mass: " << mass << endl;
866  // cerr << "Prob Density: " << density << endl;
867  prob += density;
868  }
869 
870  file.Close();
871  // cerr << "Prob: " << prob <<endl;
872  return prob;
873 }
874 
875 void listMCParticles(LCCollection *colMC) {
876  // a detail look into MC particles and their decay chain
877  if (!colMC) {
878  std::cerr << "No MC Collection Found!" << std::endl;
879  return ;
880  }
881  Int_t nMCP = colMC->getNumberOfElements();
882 
883  std::cout << setw(6) << "Index" << setw(7) << "PDG" << setw(7) << "Mother" << setw(10) << "Charge" << setw(10) << "Mass"
884  << setw(11) << "Energy" << setw(18) << "NumberOfDaughters" << setw(16) << "NumberOfParents"
885  << setw(10) << "Original" << setw(8) << "Overlay" << setw(8) << "FromSim" << std::endl;
886 
887  for (Int_t i=0;i<nMCP;i++) {
888  MCParticle *mcPart = dynamic_cast<MCParticle*>(colMC->getElementAt(i));
889  Int_t pdg = mcPart->getPDG();
890  Int_t nparents = mcPart->getParents().size();
891  Int_t motherpdg = 0;
892  if (nparents > 0) {
893  MCParticle *mother = mcPart->getParents()[0];
894  motherpdg = mother->getPDG();
895  }
896  Double_t charge = mcPart->getCharge();
897  Double_t mass = mcPart->getMass();
898  Double_t energy = mcPart->getEnergy();
899  // TVector3 pv = TVector3(mcPart->getMomentum());
900  Int_t ndaughters = mcPart->getDaughters().size();
901 #if 0
902  Int_t daughterpdg = 0;
903  if (ndaughters > 0) {
904  MCParticle *daughter = mcPart->getDaughters()[0];
905  daughterpdg = daughter->getPDG();
906  }
907 #endif
908  // TLorentzVector lortz = TLorentzVector(pv,energy);
909  Int_t originalPDG = getOriginalPDG(mcPart);
910  Int_t ioverlay = mcPart->isOverlay()? 1 : 0;
911  Int_t icreatedinsim = mcPart->isCreatedInSimulation()? 1 : 0;
912  std::cout << setw(6) << i << setw(7) << pdg << setw(7) << motherpdg << setw(10) << setprecision(3) << charge << setw(10) << mass
913  << setw(11) << energy << setw(18) << ndaughters << setw(16) << nparents
914  << setw(10) << originalPDG << setw(8) << ioverlay << setw(8) << icreatedinsim << std::endl;
915  }
916 }
917 
918  Bool_t isOverlay(ReconstructedParticle *pfo, LCCollection *colMCTL) {
919 
920  Bool_t iOverlay = false;
921 
922  LCRelationNavigator *navMCTL = new LCRelationNavigator(colMCTL);
923  LCObjectVec vecMCTL = navMCTL->getRelatedToObjects(pfo);
924  if (vecMCTL.size() > 0) {
925  MCParticle *mcPart = dynamic_cast<MCParticle *>(vecMCTL[0]);
926  if (mcPart->isOverlay()) iOverlay = true;
927  }
928 
929  return iOverlay;
930 
931 }
932 
933  void dumpJetParticles(ReconstructedParticle *jet, LCCollection *colMC, LCCollection *colMCTL) {
934  // list the particles of jets
935 
936  LCRelationNavigator *navMCTL = new LCRelationNavigator(colMCTL);
937 
938  TLorentzVector lortz = TLorentzVector(jet->getMomentum(),jet->getEnergy());
939  cerr << "----Particles of the Jet: (E,Px,Py,Pz) = (" << lortz.E() << "," << lortz.Px() << ","
940  << lortz.Py() << "," << lortz.Pz() << ")------" << endl;
941 
942  Int_t np = jet->getParticles().size();
943  cerr << setw(6) << "Index" << setw(7) << "PDG" << setw(7) << "Mother" << setw(10) << "Charge" << setw(10) << "Mass"
944  << setw(11) << "Energy" << setw(18) << "NumberOfDaughters" << setw(16) << "NumberOfParents"
945  << setw(10) << "Original" << setw(8) << "Overlay" << setw(8) << "FromSim" << endl;
946  for (Int_t i=0;i<np;i++) {
947  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle *>(jet->getParticles()[i]);
948  LCObjectVec vecMCTL = navMCTL->getRelatedToObjects(pfo);
949  if (vecMCTL.size() > 0) {
950  MCParticle *mcPart = dynamic_cast<MCParticle *>(vecMCTL[0]);
951  Int_t pdg = mcPart->getPDG();
952  Int_t nparents = mcPart->getParents().size();
953  Int_t motherpdg = 0;
954  if (nparents > 0) {
955  MCParticle *mother = mcPart->getParents()[0];
956  motherpdg = mother->getPDG();
957  }
958  Double_t charge = mcPart->getCharge();
959  Double_t mass = mcPart->getMass();
960  Double_t energy = mcPart->getEnergy();
961  TVector3 pv = TVector3(mcPart->getMomentum());
962  Int_t ndaughters = mcPart->getDaughters().size();
963 #if 0
964  Int_t daughterpdg = 0;
965  if (ndaughters > 0) {
966  MCParticle *daughter = mcPart->getDaughters()[0];
967  daughterpdg = daughter->getPDG();
968  }
969 #endif
970  Int_t originalPDG = getOriginalPDG(mcPart);
971  Int_t ioverlay = mcPart->isOverlay()? 1 : 0;
972  Int_t icreatedinsim = mcPart->isCreatedInSimulation()? 1 : 0;
973  cerr << setw(6) << i << setw(7) << pdg << setw(7) << motherpdg << setw(10) << setprecision(3) << charge
974  << setw(10) << mass << setw(11) << energy << setw(18) << ndaughters << setw(16) << nparents
975  << setw(10) << originalPDG << setw(8) << ioverlay << setw(8) << icreatedinsim << endl;
976  }
977  }
978  }
979 
980  // *******************************************************
981  // *******************************************************
982  Int_t getVertexComponents(LCCollection *colMC, LCCollection *colMCTL, ReconstructedParticle *vertex, Double_t energyComponents[2], Int_t nparticles[2]) {
983  // calculate energies from signal particles or overlay in a vertex
984  // [0]: from signal; [1]: from overlay
985 
986  for (Int_t i=0;i<2;i++) {
987  energyComponents[i] = 0.;
988  nparticles[i] = 0;
989  }
990 
991  LCRelationNavigator *navMCTL = new LCRelationNavigator(colMCTL);
992 
993  Int_t np = vertex->getParticles().size();
994  for (Int_t i=0;i<np;i++) {
995  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle *>(vertex->getParticles()[i]);
996  Double_t energy = pfo->getEnergy();
997  LCObjectVec vecMCTL = navMCTL->getRelatedToObjects(pfo);
998  if (vecMCTL.size() > 0) {
999  MCParticle *mcPart = dynamic_cast<MCParticle *>(vecMCTL[0]);
1000  if (mcPart->isOverlay()) { // overlay
1001  energyComponents[1] += energy;
1002  nparticles[1] += 1;
1003  }
1004  else {
1005  energyComponents[0] += energy;
1006  nparticles[0] += 1;
1007  }
1008  }
1009  }
1010 
1011  // find the largest energy component
1012  if (energyComponents[0] > energyComponents[1]) {
1013  return 0;
1014  }
1015  else {
1016  return 1;
1017  }
1018  }
1019 
1020  // *******************************************************
1021  // *******************************************************
1022  Int_t isVertexFromOverlay( ReconstructedParticle *vertex, LCCollection *colMC, LCCollection *colMCTL) {
1023  // 1: from overlay;
1024  // 0: from signal;
1025  Double_t energy[2];
1026  Int_t np[2];
1027  return getVertexComponents(colMC,colMCTL,vertex,energy,np);
1028 
1029  }
1030 
1031  // *******************************************************
1032  // *******************************************************
1033  Double_t getJetDistance( ReconstructedParticle *i, ReconstructedParticle *j, TString algorithm, Double_t R) {
1034  // distance betweet jet i and j
1035 
1036  TLorentzVector lortz_i = TLorentzVector(i->getMomentum(),i->getEnergy());
1037  TLorentzVector lortz_j = TLorentzVector(j->getMomentum(),j->getEnergy());
1038  Double_t pt_i = lortz_i.Pt();
1039  Double_t pt_j = lortz_j.Pt();
1040  Double_t ptmin = pt_i < pt_j ? pt_i : pt_j;
1041  Double_t y = 999.;
1042  if (algorithm == "kt") {
1043  Double_t phi_i = lortz_i.Phi();
1044  Double_t phi_j = lortz_j.Phi();
1045  Double_t rapidity_i = lortz_i.Rapidity();
1046  Double_t rapidity_j = lortz_j.Rapidity();
1047  Double_t deltaR2 = (phi_i-phi_j)*(phi_i-phi_j)+(rapidity_i-rapidity_j)*(rapidity_i-rapidity_j);
1048  y = ptmin*ptmin*deltaR2/R/R;
1049  }
1050  return y;
1051  }
1052 
1053  // *******************************************************
1054  // *******************************************************
1055  Int_t isFromVertex( ReconstructedParticle *recPart, LCCollection *colVertex) {
1056  // check if the particle is from any reconstructed vertex
1057 
1058  Int_t nVtx = colVertex->getNumberOfElements();
1059  for (Int_t i=0;i<nVtx;i++) {
1060  ReconstructedParticle *vertex = dynamic_cast<ReconstructedParticle*>(colVertex->getElementAt(i));
1061  Int_t np = vertex->getParticles().size();
1062  for (Int_t j=0;j<np;j++) {
1063  ReconstructedParticle *pfo = dynamic_cast<ReconstructedParticle *>(vertex->getParticles()[j]);
1064  if (pfo == recPart) return 1;
1065  }
1066  }
1067 
1068  return 0;
1069  }
1070 
1071  // *******************************************************
1072  // *******************************************************
1073  std::vector<Int_t> getHiggsDecayModes(LCCollection *colMC) {
1074  // get the Higgs decay mode by MC information
1075 
1076  std::vector<Int_t> nHDecay;
1077  Int_t nHbb=0,nHWW=0,nWqq=0,nWlv=0,nHgg=0,nHZZ=0,nHcc=0;
1078  Int_t nHmumu=0,nHgamgam=0,nZvv=0,nZqq=0,nZee=0,nZmm=0,nHtautau=0;
1079  Int_t nMCP = colMC->getNumberOfElements();
1080  for (Int_t i=0;i<nMCP;i++) {
1081  MCParticle *mcPart = dynamic_cast<MCParticle*>(colMC->getElementAt(i));
1082  Int_t pdg = mcPart->getPDG();
1083  Int_t ndaughters = mcPart->getDaughters().size();
1084  Int_t daughterpdg = 0;
1085  if (ndaughters > 0) {
1086  MCParticle *daughter = mcPart->getDaughters()[0];
1087  daughterpdg = daughter->getPDG();
1088  }
1089  Int_t nparents = mcPart->getParents().size();
1090  Int_t motherpdg = 0;
1091  if (nparents > 0) {
1092  MCParticle *mother = mcPart->getParents()[0];
1093  motherpdg = mother->getPDG();
1094  }
1095  if (pdg == 25 && abs(daughterpdg) == 5) nHbb++;
1096  if (pdg == 25 && abs(daughterpdg) == 24)nHWW++;
1097  if (pdg == 25 && abs(daughterpdg) == 21) nHgg++;
1098  if (pdg == 25 && abs(daughterpdg) == 23) nHZZ++;
1099  if (pdg == 25 && abs(daughterpdg) == 4) nHcc++;
1100  if (pdg == 25 && abs(daughterpdg) == 13) nHmumu++;
1101  if (pdg == 25 && abs(daughterpdg) == 15) nHtautau++;
1102  if (pdg == 25 && abs(daughterpdg) == 22) nHgamgam++;
1103  if (abs(pdg) == 24 && abs(motherpdg) == 25 && abs(daughterpdg) > 10 && abs(daughterpdg) < 20) nWlv++;
1104  if (abs(pdg) == 24 && abs(motherpdg) == 25 && abs(daughterpdg) < 10 && abs(daughterpdg) > 0) nWqq++;
1105  if (abs(pdg) == 23 && abs(motherpdg) == 25 && (abs(daughterpdg) == 12 || abs(daughterpdg) == 14 || abs(daughterpdg) == 16)) nZvv++;
1106  if (abs(pdg) == 23 && abs(motherpdg) == 25 && abs(daughterpdg) < 10 && abs(daughterpdg) > 0) nZqq++;
1107  if (abs(pdg) == 23 && abs(motherpdg) == 25 && abs(daughterpdg) == 11) nZee++;
1108  if (abs(pdg) == 23 && abs(motherpdg) == 25 && abs(daughterpdg) == 13) nZmm++;
1109  }
1110 
1111  nHDecay.push_back(nHbb);
1112  nHDecay.push_back(nHWW);
1113  nHDecay.push_back(nHgg);
1114  nHDecay.push_back(nHtautau);
1115  nHDecay.push_back(nHcc);
1116  nHDecay.push_back(nHZZ);
1117  nHDecay.push_back(nHgamgam);
1118  nHDecay.push_back(nHmumu);
1119  nHDecay.push_back(nWlv);
1120  nHDecay.push_back(nWqq);
1121  nHDecay.push_back(nZvv);
1122  nHDecay.push_back(nZqq);
1123  nHDecay.push_back(nZee);
1124  nHDecay.push_back(nZmm);
1125  return nHDecay;
1126  }
1127 
1128  // *******************************************************
1129  // *******************************************************
1130  std::vector<Int_t> getNumberOfOverlayEvents(Double_t fEcm, LCCollection *colMC) {
1131  // return number of overlay event per signal event by MC
1132 
1133  std::vector<Int_t> nOvlEvents;
1134  Int_t nOvl=0,nBgam=0,nBelectron=0,nBpositron=0;
1135  Int_t nMCP = colMC->getNumberOfElements();
1136  for (Int_t i=0;i<nMCP;i++) {
1137  MCParticle *mcPart = dynamic_cast<MCParticle*>(colMC->getElementAt(i));
1138  Int_t pdg = mcPart->getPDG();
1139 #if 0
1140  Int_t ndaughters = mcPart->getDaughters().size();
1141  Int_t daughterpdg = 0;
1142  if (ndaughters > 0) {
1143  MCParticle *daughter = mcPart->getDaughters()[0];
1144  daughterpdg = daughter->getPDG();
1145  }
1146  Int_t nparents = mcPart->getParents().size();
1147  Int_t motherpdg = 0;
1148  if (nparents > 0) {
1149  MCParticle *mother = mcPart->getParents()[0];
1150  motherpdg = mother->getPDG();
1151  }
1152 #endif
1153  Int_t ibeam = TMath::Abs(mcPart->getEnergy()-fEcm/2) < 0.1 ? 1 : 0;
1154  Int_t iovl = mcPart->isOverlay() ? 1 : 0;
1155  if (iovl == 1 && ibeam == 1 && pdg == 22) nBgam +=1;
1156  if (iovl == 1 && ibeam == 1 && pdg == 11) nBelectron +=1;
1157  if (iovl == 1 && ibeam == 1 && pdg == -11) nBpositron +=1;
1158  }
1159 
1160  nOvl = (nBgam + nBelectron + nBpositron)/2;
1161  nOvlEvents.push_back(nOvl);
1162  nOvlEvents.push_back(nBgam);
1163  nOvlEvents.push_back(nBelectron);
1164  nOvlEvents.push_back(nBpositron);
1165 
1166  return nOvlEvents;
1167  }
1168 
1169  // *******************************************************
1170  // *******************************************************
1171  void mcDebug(LCCollection *colMC) {
1172  // list the generator particles
1173 
1174  Int_t nMCP = colMC->getNumberOfElements();
1175  cerr << setw(6) << "Index" << setw(7) << "PDG" << setw(7) << "Mother" << setw(10) << "Charge" << setw(10) << "Mass" << setw(11) << "Energy" << setw(18) << "NumberOfDaughters" << setw(16) << "NumberOfParents" << setw(10) << "Original" << setw(8) << "Overlay" << setw(8) << "FromSim" << endl;
1176  for (Int_t i=0;i<nMCP;i++) {
1177  MCParticle *mcPart = dynamic_cast<MCParticle*>(colMC->getElementAt(i));
1178  Int_t pdg = mcPart->getPDG();
1179  Int_t nparents = mcPart->getParents().size();
1180  Int_t motherpdg = 0;
1181  if (nparents > 0) {
1182  MCParticle *mother = mcPart->getParents()[0];
1183  motherpdg = mother->getPDG();
1184  }
1185  Double_t charge = mcPart->getCharge();
1186  Double_t mass = mcPart->getMass();
1187  Double_t energy = mcPart->getEnergy();
1188  TVector3 pv = TVector3(mcPart->getMomentum());
1189  Int_t ndaughters = mcPart->getDaughters().size();
1190 #if 0
1191  Int_t daughterpdg = 0;
1192  if (ndaughters > 0) {
1193  MCParticle *daughter = mcPart->getDaughters()[0];
1194  daughterpdg = daughter->getPDG();
1195  }
1196 #endif
1197  TLorentzVector lortz = TLorentzVector(pv,energy);
1198  Int_t originalPDG = getOriginalPDG(mcPart);
1199  Int_t ioverlay = mcPart->isOverlay()? 1 : 0;
1200  Int_t icreatedinsim = mcPart->isCreatedInSimulation()? 1 : 0;
1201  cerr << setw(6) << i << setw(7) << pdg << setw(7) << motherpdg << setw(10) << setprecision(3) << charge << setw(10) << mass << setw(11) << energy << setw(18) << ndaughters << setw(16) << nparents << setw(10) << originalPDG << setw(8) << ioverlay << setw(8) << icreatedinsim << endl;
1202  }
1203 
1204  }
1205 
1206  // *******************************************************
1207  // *******************************************************
1208  TLorentzVector getLorentzEcm(Double_t fEcm) {
1209  const Double_t fCrossAngle = 0.014;
1210  return TLorentzVector(fEcm*TMath::Sin(fCrossAngle/2),0.,0.,fEcm);
1211  }
1212  TLorentzVector getLorentzEcm(Double_t fEcm,Bool_t isCrossingAngle) {
1213  const Double_t fCrossAngle = 0.014;
1214  if (isCrossingAngle) {
1215  return TLorentzVector(fEcm*TMath::Sin(fCrossAngle/2),0.,0.,fEcm);
1216  }
1217  else {
1218  return TLorentzVector(0.,0.,0.,fEcm);
1219  }
1220  }
1221 
1222  // *******************************************************
1223  // *******************************************************
1224  Double_t getRecoilMass(TLorentzVector lortzEcm, TLorentzVector lortzZ) {
1225  TLorentzVector lortzRecoil = lortzEcm - lortzZ;
1226  return lortzRecoil.M();
1227  }
1228 
1229  // *******************************************************
1230  // *******************************************************
1231  Double_t getAcoPlanarity(TLorentzVector lortz1, TLorentzVector lortz2) {
1232  Double_t phi1 = lortz1.Phi();
1233  Double_t phi2 = lortz2.Phi();
1234  Double_t acop12 = TMath::Abs(phi1-phi2);
1235  return acop12 > TMath::Pi() ? TMath::Pi()*2. - acop12 : acop12;
1236  }
1237 
1238  // *******************************************************
1239  // *******************************************************
1240  TLorentzVector getLorentzVector(ReconstructedParticle *pfo) {
1241  return TLorentzVector(pfo->getMomentum(),pfo->getEnergy());
1242  }
1243 
1244  // *******************************************************
1245  // *******************************************************
1246  Double_t getCosTheta(ReconstructedParticle *part1, ReconstructedParticle *part2) {
1247  TVector3 momentum1 = TVector3(part1->getMomentum());
1248  TVector3 momentum2 = TVector3(part2->getMomentum());
1249  Double_t cosTheta = momentum1.Dot(momentum2)/momentum1.Mag()/momentum2.Mag();
1250  return cosTheta;
1251  }
1252 
1253  // *******************************************************
1254  // *******************************************************
1255  Double_t getCosTheta(TLorentzVector part1, TLorentzVector part2) {
1256  TVector3 momentum1 = TVector3(part1.Vect());
1257  TVector3 momentum2 = TVector3(part2.Vect());
1258  Double_t cosTheta = momentum1.Dot(momentum2)/momentum1.Mag()/momentum2.Mag();
1259  return cosTheta;
1260  }
1261 
1262 } // namespace mylib
isolep::getLeptonID
Int_t getLeptonID(ReconstructedParticle *recPart)
Definition: Utilities.cc:237
isolep::getSplitTag
Bool_t getSplitTag(ReconstructedParticle *motherPart, ReconstructedParticle *recPart)
Definition: Utilities.cc:308
isolep::getOriginalPDGForIsoLep
Int_t getOriginalPDGForIsoLep(MCParticle *mcPart, LCCollection *colMC)
Definition: Utilities.cc:150
isolep::mcDebug
void mcDebug(LCCollection *colMC)
Definition: Utilities.cc:1171
isolep::getLikelihood
Double_t getLikelihood(TString fname, TString hist, Double_t mass)
Definition: Utilities.cc:849
isolep::getOriginalSerialForZHH
Int_t getOriginalSerialForZHH(MCParticle *mcPart, LCCollection *colMCP)
Definition: Utilities.cc:214
isolep::getCosTheta
Double_t getCosTheta(ReconstructedParticle *part1, ReconstructedParticle *part2)
Definition: Utilities.cc:1246
isolep::getJetDistance
Double_t getJetDistance(ReconstructedParticle *i, ReconstructedParticle *j, TString algorithm, Double_t R)
Definition: Utilities.cc:1033
isolep::isVertexFromOverlay
Int_t isVertexFromOverlay(ReconstructedParticle *vertex, LCCollection *colMC, LCCollection *colMCTL)
Definition: Utilities.cc:1022
isolep::calculateEnergyComponents
Int_t calculateEnergyComponents(LCCollection *colMC, LCCollection *colMCTL, ReconstructedParticle *jet, Double_t energyComponents[4])
Definition: Utilities.cc:745
isolep::dumpJetParticles
void dumpJetParticles(ReconstructedParticle *jet, LCCollection *colMC, LCCollection *colMCTL)
Definition: Utilities.cc:933
isolep::getConeEnergy
Double_t getConeEnergy(ReconstructedParticle *recPart, LCCollection *colPFO, Double_t cosCone)
Definition: Utilities.cc:361
isolep::getFSRMomentum
TLorentzVector getFSRMomentum(ReconstructedParticle *recPart, LCCollection *colPFO)
Definition: Utilities.cc:469
isolep::getLinkedMCParticle
Int_t getLinkedMCParticle(ReconstructedParticle *recPart, LCCollection *colMCTL, Double_t &weight, Int_t &nMCTL)
Definition: Utilities.cc:61
isolep::getNumberOfOverlayEvents
std::vector< Int_t > getNumberOfOverlayEvents(Double_t fEcm, LCCollection *colMC)
Definition: Utilities.cc:1130
isolep::getRecoilMass
Double_t getRecoilMass(TLorentzVector lortzEcm, TLorentzVector lortzZ)
Definition: Utilities.cc:1224
isolep::getAcoPlanarity
Double_t getAcoPlanarity(TLorentzVector lortz1, TLorentzVector lortz2)
Definition: Utilities.cc:1231
isolep::doPhotonRecovery
void doPhotonRecovery(ReconstructedParticle *electron, LCCollection *colPFO, ReconstructedParticleImpl *recoElectron, Double_t fCosFSRCut=0.999)
Definition: Utilities.cc:645
isolep::getOriginalPDG
Int_t getOriginalPDG(MCParticle *mcPart, Bool_t iHiggs=0)
Definition: Utilities.cc:85
isolep::isFromVertex
Int_t isFromVertex(ReconstructedParticle *recPart, LCCollection *colVertex)
Definition: Utilities.cc:1055
isolep::getLorentzEcm
TLorentzVector getLorentzEcm(Double_t fEcm)
Definition: Utilities.cc:1208
isolep::getMCSerial
Int_t getMCSerial(MCParticle *mcPart, LCCollection *colMCP)
Definition: Utilities.cc:16
isolep::getInvariantMass
Double_t getInvariantMass(ReconstructedParticle *recPart1, ReconstructedParticle *recPart2)
Definition: Utilities.cc:577
isolep::getLorentzVector
TLorentzVector getLorentzVector(ReconstructedParticle *pfo)
Definition: Utilities.cc:1240
isolep::getMCParticle
MCParticle * getMCParticle(ReconstructedParticle *recPart, LCCollection *colMCTL)
Definition: Utilities.cc:30
isolep::isOverlay
Bool_t isOverlay(ReconstructedParticle *pfo, LCCollection *colMCTL)
Definition: Utilities.cc:918
isolep::getVertexComponents
Int_t getVertexComponents(LCCollection *colMC, LCCollection *colMCTL, ReconstructedParticle *vetex, Double_t energyComponents[2], Int_t nparticles[2])
Definition: Utilities.cc:982
isolep::getHiggsDecayModes
std::vector< Int_t > getHiggsDecayModes(LCCollection *colMC)
Definition: Utilities.cc:1073
isolep::getOriginalSerial
Int_t getOriginalSerial(MCParticle *mcPart, LCCollection *colMCP, Bool_t iHiggs=0)
Definition: Utilities.cc:179
isolep::getFSRTag
Bool_t getFSRTag(ReconstructedParticle *motherPart, ReconstructedParticle *recPart, Double_t fCosFSRCut=0.999)
Definition: Utilities.cc:289
isolep::isFoundInVector
Bool_t isFoundInVector(ReconstructedParticle *pfo, std::vector< lcio::ReconstructedParticle * > &pfos)
Definition: Utilities.cc:706
isolep::jetFunction
Double_t jetFunction(TLorentzVector lortz, float beta, float power)
Definition: Utilities.cc:715
isolep::isSelectedByFastJet
Int_t isSelectedByFastJet(ReconstructedParticle *pfo, LCCollection *colFastJet, Double_t &ratioEPartEJet, Double_t &ratioPTMJet)
Definition: Utilities.cc:613
isolep::listMCParticles
void listMCParticles(LCCollection *colMC)
Definition: Utilities.cc:875