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.7 by ceballos, Tue Jan 18 16:41:13 2011 UTC

# Line 15 | Line 15 | JetTools::~JetTools()
15   }
16  
17   //Remember to remove the signal from particles before inputting into the function
18 < Double_t JetTools::NJettiness(const ParticleOArr *particles, const JetOArr *jets, bool UseQ, Double_t Y){
18 > Double_t JetTools::NJettiness(const ParticleOArr *particles, const JetOArr *jets, double Q, double Y){
19    if(particles->GetEntries() <= 0) return 0.0;
20  
21    Double_t fval = 0.0;
# Line 32 | Line 32 | Double_t JetTools::NJettiness(const Part
32      fval = fval + fvalpart;
33    }
34  
35 <  if(UseQ == kTRUE) fval = fval / particles->At(0)->Pt();
35 >  fval = fval / Q;
36  
37    return fval;
38   }
39  
40 < Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, bool UseQ, Double_t Y){
40 > Double_t JetTools::NJettiness(const PFCandidateOArr *pfCandidates, const JetOArr *jets, double Q, double Y){
41 >  if(pfCandidates->GetEntries() <= 0) return 0.0;
42 >
43 >  Double_t fval = 0.0;
44 >  Double_t fvalpart;
45 >
46 >  for(int i=0;i<int(pfCandidates->GetEntries());i++){
47 >    fvalpart = (pfCandidates->At(i)->Pt()) * TMath::Exp(-TMath::Abs(pfCandidates->At(i)->Eta()-Y));
48 >
49 >    for(int j=0;j<int(jets->GetEntries());j++){
50 >      fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
51 >                 (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-pfCandidates->At(i)->Eta()))
52 >                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi()))));
53 >    }
54 >    fval = fval + fvalpart;
55 >  }
56 >
57 >  fval = fval / Q;
58 >
59 >  return fval;
60 > }
61 >
62 > Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, double Q, double Y){
63    if(tracks->GetEntries() <= 0) return 0.0;
64  
65    Double_t fval = 0.0;
# Line 54 | Line 76 | Double_t JetTools::NJettiness(const Trac
76      fval = fval + fvalpart;
77    }
78  
79 <  if(UseQ == kTRUE) fval = fval / tracks->At(0)->Pt();
79 >  fval = fval / Q;
80    
81    return fval;
82   }
83  
84 < Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, bool UseQ, Double_t Y){
84 > Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, double Q, double Y){
85    if(jetsS->GetEntries() <= 0) return 0.0;
86  
87    Double_t fval = 0.0;
# Line 76 | Line 98 | Double_t JetTools::NJettiness(const JetO
98      fval = fval + fvalpart;
99    }
100  
101 <  if(UseQ == kTRUE) fval = fval / jetsS->At(0)->Pt();
101 >  fval = fval / Q;
102    
103    return fval;
104   }
105  
106 < Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, bool UseQ, Double_t Y){
106 > Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, double Q, double Y){
107    if(calos->GetEntries() <= 0) return 0.0;
108  
109    Double_t fval = 0.0;
# Line 98 | Line 120 | Double_t JetTools::NJettiness(const Calo
120      fval = fval + fvalpart;
121    }
122  
123 <  if(UseQ == kTRUE) fval = fval / calos->At(0)->Pt();
123 >  fval = fval / Q;
124    
125    return fval;
126   }
# Line 184 | Line 206 | Double_t JetTools::CosineOmega(const Par
206  
207   //Transverse Higgs mass
208   Double_t JetTools::MtHiggs(const ParticleOArr * leptons,
209 <                           const Met *met, int nsel){
209 >                           const Met *met, double metFraction[2], int nsel){
210    if(leptons->Entries() < 2) return -999.0;
211  
212    double mtHiggs = -999.0;
# Line 230 | Line 252 | Double_t JetTools::MtHiggs(const Particl
252      mll   = dilepton->Mass();
253      mnu   = 0.0;
254    }
255 <  else if(nsel == 4){ // Use the formula from hep-ph:1006.4998
255 >  else if(nsel == 4){ // Use of Mt mass and replacing mnu using the met optimal
256 >    enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
257 >    enenn = TMath::Sqrt(met->Pt() *met->Pt()  + 0.0*0.0);
258 >    enex  = dilepton->Px() + met->Px();
259 >    eney  = dilepton->Py() + met->Py();
260 >    mll   = dilepton->Mass();
261 >    double metAuxPx[2] = {met->Px() * metFraction[0],
262 >                          met->Px() * (1.0 - metFraction[0])};
263 >    double metAuxPy[2] = {met->Py() * metFraction[1],
264 >                          met->Py() * (1.0 - metFraction[1])};
265 >    double ene = TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
266 >                 TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
267 >    double px = metAuxPx[0] + metAuxPx[1];
268 >    double py = metAuxPy[0] + metAuxPy[1];
269 >    mnu = TMath::Sqrt(ene*ene - px*px - py*py);
270 >  }
271 >  else if(nsel == 5){ // Using the optimal met value
272 >    double metAuxPx[2] = {met->Px() * metFraction[0],
273 >                          met->Px() * (1.0 - metFraction[0])};
274 >    double metAuxPy[2] = {met->Py() * metFraction[1],
275 >                          met->Py() * (1.0 - metFraction[1])};
276 >    double ene = leptons->At(0)->Pt() + leptons->At(1)->Pt() +
277 >                 TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
278 >                 TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
279 >    double px = leptons->At(0)->Px() + leptons->At(1)->Px() +
280 >                metAuxPx[0] + metAuxPx[1];
281 >    double py = leptons->At(0)->Py() + leptons->At(1)->Py() +
282 >                metAuxPy[0] + metAuxPy[1];
283 >    mtHiggs = ene*ene - px*px - py*py;
284 >  }
285 >  else if(nsel == 6){ // Use the formula from hep-ph:1006.4998
286      mtHiggs = 2*leptons->At(0)->Pt()*leptons->At(0)->Pt() + 2*leptons->At(1)->Pt()*leptons->At(1)->Pt() + 3 * (
287        leptons->At(0)->Pt()*leptons->At(1)->Pt() + met->Pt()*(leptons->At(0)->Pt()+leptons->At(1)->Pt())
288        - met->Px()*dilepton->Px() - met->Py()*dilepton->Py()
289        - leptons->At(0)->Px()*leptons->At(1)->Px() - leptons->At(0)->Py()*leptons->At(1)->Py());
290    }
291  
292 <  if(nsel >= 0 && nsel <= 3){
292 >  if(nsel >= 0 && nsel <= 4){
293      mtHiggs = mll*mll + mnu*mnu + 2.0*(enell*enenn - enex*enex - eney*eney);
294    }
295    if(mtHiggs <= 0) mtHiggs = 0.0;
# Line 247 | Line 299 | Double_t JetTools::MtHiggs(const Particl
299  
300    return mtHiggs;
301   }
302 +
303 + void JetTools::Alpha(Double_t AlphaVar[2], const TrackCol *tracks, Jet *jet, const VertexCol *vertices, Double_t  delta_z, Double_t delta_cone){  
304 +  AlphaVar[0] = -1.0;
305 +  AlphaVar[1] = -1.0;
306 +  if(tracks->GetEntries() <= 0) return;
307 +
308 +  double Pt_jets_X = 0. ;
309 +  double Pt_jets_Y = 0. ;
310 +  double Pt_jets_X_tot = 0. ;
311 +  double Pt_jets_Y_tot = 0. ;
312 +
313 +  for(int i=0;i<int(tracks->GetEntries());i++){
314 +    if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
315 +      Pt_jets_X_tot += tracks->At(i)->Px();
316 +      Pt_jets_Y_tot += tracks->At(i)->Py();  
317 +      double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertices->At(0)));
318 +      if(pDz < delta_z){
319 +        Pt_jets_X += tracks->At(i)->Px();
320 +        Pt_jets_Y += tracks->At(i)->Py();  
321 +      }
322 +    }
323 +  }
324 +
325 +  if(jet->Pt() > 0)
326 +    AlphaVar[0] = sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / jet->Pt();
327 +  if(Pt_jets_X_tot > 0 || Pt_jets_Y_tot > 0)
328 +    AlphaVar[1] = 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);
329 + }
330 +
331 +
332 + void JetTools::Alpha(Double_t AlphaVar[2], const PFJet *jet, const VertexCol *vertices, Double_t  delta_z){  
333 +  AlphaVar[0] = -1.0;
334 +  AlphaVar[1] = -1.0;
335 +
336 +  double Pt_jets_X = 0. ;
337 +  double Pt_jets_Y = 0. ;
338 +  double Pt_jets_X_tot = 0. ;
339 +  double Pt_jets_Y_tot = 0. ;
340 +
341 +  for(UInt_t i=0;i<jet->NPFCands();i++){
342 +    if(jet->PFCand(i)->BestTrk()){
343 +      Pt_jets_X_tot += jet->PFCand(i)->BestTrk()->Px();
344 +      Pt_jets_Y_tot += jet->PFCand(i)->BestTrk()->Py();  
345 +      double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(0)));
346 +      if(pDz < delta_z){
347 +        Pt_jets_X += jet->PFCand(i)->BestTrk()->Px();
348 +        Pt_jets_Y += jet->PFCand(i)->BestTrk()->Py();  
349 +      }
350 +    }
351 +  }
352 +
353 +  if(jet->Pt() > 0)
354 +    AlphaVar[0] = sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / jet->Pt();
355 +  if(Pt_jets_X_tot > 0 || Pt_jets_Y_tot > 0)
356 +    AlphaVar[1] = 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);
357 + }
358 +

Diff Legend

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