ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetTools.cc
(Generate patch)

Comparing UserCode/MitPhysics/Utils/src/JetTools.cc (file contents):
Revision 1.3 by ceballos, Fri Jul 30 20:23:31 2010 UTC vs.
Revision 1.14 by ceballos, Wed Apr 20 10:06:15 2011 UTC

# Line 1 | Line 1
1   #include "MitPhysics/Utils/interface/JetTools.h"
2 + #include <algorithm>
3 + #include <vector>
4  
5   ClassImp(mithep::JetTools)
6  
7   using namespace mithep;
8 +
9  
10   JetTools::JetTools()
11   {
# Line 15 | Line 18 | JetTools::~JetTools()
18   }
19  
20   //Remember to remove the signal from particles before inputting into the function
21 < Double_t JetTools::NJettiness(const ParticleOArr *particles, const JetOArr *jets, bool UseQ, Double_t Y){
21 > Double_t JetTools::NJettiness(const ParticleOArr *particles, const JetOArr *jets, double Q, double Y){
22    if(particles->GetEntries() <= 0) return 0.0;
23  
24    Double_t fval = 0.0;
# Line 32 | Line 35 | Double_t JetTools::NJettiness(const Part
35      fval = fval + fvalpart;
36    }
37  
38 <  if(UseQ == kTRUE) fval = fval / particles->At(0)->Pt();
38 >  fval = fval / Q;
39  
40    return fval;
41   }
42  
43 < Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, bool UseQ, Double_t Y){
43 > Double_t JetTools::NJettiness(const PFCandidateOArr *pfCandidates, const JetOArr *jets, double Q, double Y){
44 >  if(pfCandidates->GetEntries() <= 0) return 0.0;
45 >
46 >  Double_t fval = 0.0;
47 >  Double_t fvalpart;
48 >
49 >  for(int i=0;i<int(pfCandidates->GetEntries());i++){
50 >    fvalpart = (pfCandidates->At(i)->Pt()) * TMath::Exp(-TMath::Abs(pfCandidates->At(i)->Eta()-Y));
51 >
52 >    for(int j=0;j<int(jets->GetEntries());j++){
53 >      fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
54 >                 (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-pfCandidates->At(i)->Eta()))
55 >                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi()))));
56 >    }
57 >    fval = fval + fvalpart;
58 >  }
59 >
60 >  fval = fval / Q;
61 >
62 >  return fval;
63 > }
64 >
65 > Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, double Q, double Y){
66    if(tracks->GetEntries() <= 0) return 0.0;
67  
68    Double_t fval = 0.0;
# Line 54 | Line 79 | Double_t JetTools::NJettiness(const Trac
79      fval = fval + fvalpart;
80    }
81  
82 <  if(UseQ == kTRUE) fval = fval / tracks->At(0)->Pt();
82 >  fval = fval / Q;
83    
84    return fval;
85   }
86  
87 < Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, bool UseQ, Double_t Y){
87 > Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, double Q, double Y){
88    if(jetsS->GetEntries() <= 0) return 0.0;
89  
90    Double_t fval = 0.0;
# Line 76 | Line 101 | Double_t JetTools::NJettiness(const JetO
101      fval = fval + fvalpart;
102    }
103  
104 <  if(UseQ == kTRUE) fval = fval / jetsS->At(0)->Pt();
104 >  fval = fval / Q;
105    
106    return fval;
107   }
108  
109 < Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, bool UseQ, Double_t Y){
109 > Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, double Q, double Y){
110    if(calos->GetEntries() <= 0) return 0.0;
111  
112    Double_t fval = 0.0;
# Line 98 | Line 123 | Double_t JetTools::NJettiness(const Calo
123      fval = fval + fvalpart;
124    }
125  
126 <  if(UseQ == kTRUE) fval = fval / calos->At(0)->Pt();
126 >  fval = fval / Q;
127    
128    return fval;
129   }
# Line 184 | Line 209 | Double_t JetTools::CosineOmega(const Par
209  
210   //Transverse Higgs mass
211   Double_t JetTools::MtHiggs(const ParticleOArr * leptons,
212 <                           const Met *met, int nsel){
212 >                           const Met *met, double metFraction[2], int nsel){
213    if(leptons->Entries() < 2) return -999.0;
214  
215    double mtHiggs = -999.0;
# Line 230 | Line 255 | Double_t JetTools::MtHiggs(const Particl
255      mll   = dilepton->Mass();
256      mnu   = 0.0;
257    }
258 <  else if(nsel == 4){ // Use the formula from hep-ph:1006.4998
258 >  else if(nsel == 4){ // Use of Mt mass and replacing mnu using the met optimal
259 >    enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
260 >    enenn = TMath::Sqrt(met->Pt() *met->Pt()  + 0.0*0.0);
261 >    enex  = dilepton->Px() + met->Px();
262 >    eney  = dilepton->Py() + met->Py();
263 >    mll   = dilepton->Mass();
264 >    double metAuxPx[2] = {met->Px() * metFraction[0],
265 >                          met->Px() * (1.0 - metFraction[0])};
266 >    double metAuxPy[2] = {met->Py() * metFraction[1],
267 >                          met->Py() * (1.0 - metFraction[1])};
268 >    double ene = TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
269 >                 TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
270 >    double px = metAuxPx[0] + metAuxPx[1];
271 >    double py = metAuxPy[0] + metAuxPy[1];
272 >    mnu = TMath::Sqrt(ene*ene - px*px - py*py);
273 >  }
274 >  else if(nsel == 5){ // Using the optimal met value
275 >    double metAuxPx[2] = {met->Px() * metFraction[0],
276 >                          met->Px() * (1.0 - metFraction[0])};
277 >    double metAuxPy[2] = {met->Py() * metFraction[1],
278 >                          met->Py() * (1.0 - metFraction[1])};
279 >    double ene = leptons->At(0)->Pt() + leptons->At(1)->Pt() +
280 >                 TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
281 >                 TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
282 >    double px = leptons->At(0)->Px() + leptons->At(1)->Px() +
283 >                metAuxPx[0] + metAuxPx[1];
284 >    double py = leptons->At(0)->Py() + leptons->At(1)->Py() +
285 >                metAuxPy[0] + metAuxPy[1];
286 >    mtHiggs = ene*ene - px*px - py*py;
287 >  }
288 >  else if(nsel == 6){ // Use the formula from hep-ph:1006.4998
289      mtHiggs = 2*leptons->At(0)->Pt()*leptons->At(0)->Pt() + 2*leptons->At(1)->Pt()*leptons->At(1)->Pt() + 3 * (
290        leptons->At(0)->Pt()*leptons->At(1)->Pt() + met->Pt()*(leptons->At(0)->Pt()+leptons->At(1)->Pt())
291        - met->Px()*dilepton->Px() - met->Py()*dilepton->Py()
292        - leptons->At(0)->Px()*leptons->At(1)->Px() - leptons->At(0)->Py()*leptons->At(1)->Py());
293    }
294 +  else if(nsel == 7){ // Use of M mass and mnu == 0
295 +    double deltaPhiDileptonMet = MathUtils::DeltaPhi(dilepton->Phi(),
296 +                                                     met->Phi());
297 +    mtHiggs = 2.0*dilepton->Pt()*met->Pt()*(1.0 - cos(deltaPhiDileptonMet));
298 +  }
299  
300 <  if(nsel >= 0 && nsel <= 3){
300 >  if(nsel >= 0 && nsel <= 4){
301      mtHiggs = mll*mll + mnu*mnu + 2.0*(enell*enenn - enex*enex - eney*eney);
302    }
303 +
304    if(mtHiggs <= 0) mtHiggs = 0.0;
305    else             mtHiggs = TMath::Sqrt(mtHiggs);
306  
# Line 247 | Line 308 | Double_t JetTools::MtHiggs(const Particl
308  
309    return mtHiggs;
310   }
311 +
312 + Double_t JetTools::Beta(const TrackCol *tracks, Jet *jet, const Vertex *vertex, Double_t  delta_z, Double_t delta_cone){  
313 +
314 +  if(tracks->GetEntries() <= 0) return 1.0;
315 +
316 +  double Pt_jets_X = 0. ;
317 +  double Pt_jets_Y = 0. ;
318 +  double Pt_jets_X_tot = 0. ;
319 +  double Pt_jets_Y_tot = 0. ;
320 +
321 +  for(int i=0;i<int(tracks->GetEntries());i++){
322 +    if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
323 +      Pt_jets_X_tot += tracks->At(i)->Px();
324 +      Pt_jets_Y_tot += tracks->At(i)->Py();  
325 +      double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertex));
326 +      if(pDz < delta_z){
327 +        Pt_jets_X += tracks->At(i)->Px();
328 +        Pt_jets_Y += tracks->At(i)->Py();  
329 +      }
330 +    }
331 +  }
332 +
333 +  if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
334 +    return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);
335 +
336 +  return 1.0;
337 + }
338 +
339 +
340 + Double_t JetTools::Beta(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){  
341 +  double Pt_jets= 0. ;
342 +  double Pt_jetsTot = 0. ;
343 +
344 +  for(UInt_t i=0;i<jet->NPFCands();i++){
345 +    if(jet->PFCand(i)->BestTrk()){
346 +      Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt();
347 +      double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
348 +      if(pDz < delta_z){
349 +        Pt_jets += jet->PFCand(i)->BestTrk()->Pt();
350 +      }
351 +    }
352 +  }
353 +
354 +  Double_t beta = 1.0;
355 +  if (Pt_jetsTot > 0)
356 +    beta = Pt_jets/Pt_jetsTot;
357 +
358 +  return beta;
359 + }
360 +
361 + Double_t JetTools::Beta2(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){  
362 +  double Pt_jets= 0. ;
363 +  double Pt_jetsTot = 0. ;
364 +
365 +  for(UInt_t i=0;i<jet->NPFCands();i++){
366 +    if(jet->PFCand(i)->BestTrk()){
367 +      Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
368 +      double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
369 +      if(pDz < delta_z){
370 +        Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
371 +      }
372 +    }
373 +  }
374 +
375 +  Double_t beta = 1.0;
376 +  if (Pt_jetsTot > 0)
377 +    beta = Pt_jets/Pt_jetsTot;
378 +  return beta;
379 + }
380 +
381 +
382 + Bool_t  JetTools::PassBetaVertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
383 +
384 +  Bool_t passBetaCut = kTRUE;
385 +  if(vertices->GetEntries() > 0) {
386 +    Double_t Beta = JetTools::Beta(jet, referenceVertex, 0.2);
387 +    Double_t Beta_other = 0.0;
388 +    for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
389 +      if (referenceVertex == vertices->At(nv)) continue;
390 +      Double_t BetaAux = JetTools::Beta(jet, vertices->At(nv), 0.2);
391 +      if(BetaAux > Beta_other) Beta_other = BetaAux;
392 +    }
393 +    if(Beta_other > Beta) passBetaCut = kFALSE;
394 +  }
395 +
396 +  return passBetaCut;
397 +
398 + }
399 +
400 + Bool_t  JetTools::PassBeta2VertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
401 +
402 +  Bool_t passBetaCut = kTRUE;
403 +  if(vertices->GetEntries() > 0) {
404 +    Double_t Beta = JetTools::Beta2(jet, referenceVertex, 0.2);
405 +    Double_t Beta_other = 0.0;
406 +    for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
407 +      if (referenceVertex == vertices->At(nv)) continue;
408 +      Double_t BetaAux = JetTools::Beta2(jet, vertices->At(nv), 0.2);
409 +      if(BetaAux > Beta_other) Beta_other = BetaAux;
410 +    }
411 +    if(Beta_other > Beta) passBetaCut = kFALSE;
412 +  }
413 +
414 +  return passBetaCut;
415 +
416 + }
417 +
418 +
419 + Int_t JetTools::MaxBetaVertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){  
420 +  
421 +  Int_t vertexIndex = -1;
422 +  double beta = -0.1;
423 +  for (UInt_t v=0; v < vertices->GetEntries(); v++){
424 +    Double_t betaTmp = JetTools::Beta(jet, vertices->At(v), delta_z);
425 +    if (betaTmp > beta) {
426 +      beta = betaTmp;
427 +      vertexIndex = v;
428 +    }
429 +  }
430 +  return vertexIndex;
431 +
432 + }
433 +
434 + Int_t JetTools::MaxBeta2VertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){  
435 +  
436 +  Int_t vertexIndex = -1;
437 +  double beta = -0.1;
438 +  for (UInt_t v=0; v < vertices->GetEntries(); v++){
439 +    Double_t betaTmp = JetTools::Beta2(jet, vertices->At(v), delta_z);
440 +    if (betaTmp > beta) {
441 +      beta = betaTmp;
442 +      vertexIndex = v;
443 +    }
444 +  }
445 +  return vertexIndex;
446 +
447 + }
448 +
449 +
450 + Int_t JetTools::JetToPVAssociation(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){  
451 +
452 +  std::vector<float> verticesPt2(vertices->GetEntries());
453 +  for(UInt_t i=0;i<jet->NPFCands();i++){
454 +    if(jet->PFCand(i)->BestTrk()){
455 +      double minDZ = delta_z;
456 +      int trackVertexIndex = -1;
457 +      for (UInt_t v=0; v < vertices->GetEntries(); v++){
458 +        if (minDZ > TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)))) {
459 +          minDZ = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)));
460 +          trackVertexIndex = v;
461 +        }
462 +      }
463 +      if (trackVertexIndex < 0) continue;
464 +      verticesPt2[trackVertexIndex]+= jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
465 +    }
466 +  }
467 +
468 +  Int_t vertexIndex = 0;
469 +  float pt2Max = 0;
470 +  for (uint i=0; i < verticesPt2.size(); ++i){
471 +    if (pt2Max < verticesPt2[i]) {
472 +      pt2Max = verticesPt2[i];
473 +      vertexIndex = i;
474 +    }
475 +  }
476 +  return vertexIndex;
477 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines