LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
MinMassesCombJigsaw.cc
Go to the documentation of this file.
1 // RestFrames: particle physics event analysis library
3 // --------------------------------------------------------------------
4 // Copyright (c) 2014-2016, Christopher Rogan
14 // This file is part of RestFrames.
15 //
16 // RestFrames is free software; you can redistribute it and/or modify
17 // it under the terms of the GNU General Public License as published by
18 // the Free Software Foundation; either version 2 of the License, or
19 // (at your option) any later version.
20 //
21 // RestFrames is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with RestFrames. If not, see <http://www.gnu.org/licenses/>.
29 
30 #include "RestFrames/RestFrame.hh"
32 
33 namespace RestFrames {
34 
36  // MinMassesCombJigsaw class methods
38  MinMassesCombJigsaw::MinMassesCombJigsaw(const std::string& sname,
39  const std::string& stitle) :
40  CombinatoricJigsaw(sname, stitle, 2, 2) {}
41 
42  MinMassesCombJigsaw::MinMassesCombJigsaw() : CombinatoricJigsaw() {}
43 
44  MinMassesCombJigsaw::~MinMassesCombJigsaw() {}
45 
46  MinMassesCombJigsaw& MinMassesCombJigsaw::Empty(){
47  return MinMassesCombJigsaw::m_Empty;
48  }
49 
50  void MinMassesCombJigsaw::AddFrame(const RestFrame& frame, int i){
51  if(!frame) return;
52  if(!GetGroup()) return;
53 
54  ConstRestFrameList frames =
55  frame.GetListVisibleFrames()+
56  frame.GetListInvisibleFrames();
57  int N = frames.GetN();
58  for(int f = 0; f < N; f++){
59  if(GetGroup().ContainsFrame(frames[f]))
60  AddChildFrame(frames[f], i);
61  AddDependancyFrame(frames[f], i);
62  }
63  }
64 
65  void MinMassesCombJigsaw::AddFrames(const ConstRestFrameList& frames, int i){
66  int N = frames.GetN();
67  for(int f = 0; f < N; f++)
68  AddFrame(frames[f],i);
69  }
70 
71  bool MinMassesCombJigsaw::AnalyzeEvent(){
72  if(!IsSoundMind() || !GetGroup())
73  return SetSpirit(false);
74 
75  if(!InitializeCombinatoric()){
76  m_Log << LogWarning;
77  m_Log << "Problem initializing event info" << LogEnd;
78  return SetSpirit(false);
79  }
80 
81  int Ninput = GetNInputStates();
82 
83  bool DO_N3 = (Ninput >= 2) &&
84  GetDependancyFrames(0) == GetChildFrames(0) &&
85  GetDependancyFrames(1) == GetChildFrames(1);
86  if(DO_N3){
87  DO_N3 =
88  (GetNinputForChild(0) <= 1) &&
89  (GetNinputForChild(1) <= 1) &&
90  !IsNinputExclForChild(0) &&
91  !IsNinputExclForChild(1) &&
92  !IsChargeSetForChild(0) &&
93  !IsChargeSetForChild(1) &&
94  !IsChargeSetForObject(0) &&
95  !IsChargeSetForObject(1);
96  }
97  if(!DO_N3){
98  if(!CombinatoricJigsaw::LoopCombinatoric()){
99  m_Log << LogWarning;
100  m_Log << "Problem looping over combinatorics" << LogEnd;
101  return SetSpirit(false);
102  }
103  return SetSpirit(true);
104  }
105 
106  // DO N^3 calculation
107  std::vector<TLorentzVector> inputs;
108  for(int i = 0; i < Ninput; i++){
109  inputs.push_back(GetInputState(i).GetFourVector());
110  if(inputs[i].M() < 0.) inputs[i].SetVectM(inputs[i].Vect(),0.);
111  }
112 
113  // boost input vectors to CM frame
114  TLorentzVector TOT(0.,0.,0.,0.);
115  for(int i = 0; i < Ninput; i++) TOT += inputs[i];
116  TVector3 boost = TOT.BoostVector();
117  if(boost.Mag() >= 1.)
118  boost = (1.-1.e-8)/boost.Mag()*boost;
119  for(int i = 0; i < Ninput; i++) inputs[i].Boost(-boost);
120 
121  int ip_max[2];
122  int jp_max[2];
123  for(int i = 0; i < 2; i++) ip_max[i] = -1;
124  for(int i = 0; i < 2; i++) jp_max[i] = -1;
125  double metric_max = -1.;
126  // Loop over all 2-jet seed probes
127  int ip[2], jp[2];
128  for(ip[0] = 0; ip[0] < Ninput-1; ip[0]++){
129  for(ip[1] = ip[0]+1; ip[1] < Ninput; ip[1]++){
130  TVector3 nRef = inputs[ip[0]].Vect().Cross(inputs[ip[1]].Vect());
131  int Nhem[2];
132  TLorentzVector hem[2];
133  for(int i = 0; i < 2; i++){
134  Nhem[i] = 0;
135  hem[i].SetPxPyPzE(0.,0.,0.,0.);
136  }
137  // Loop over all jets
138  for(int i = 0; i < Ninput; i++){
139  if((i == ip[0]) || (i ==ip[1])) continue;
140  int ihem = int(inputs[i].Vect().Dot(nRef) > 0.);
141  Nhem[ihem]++;
142  hem[ihem] += inputs[i];
143  }
144  // assign 2 probes
145  for(jp[0] = 0; jp[0] < 2; jp[0]++){
146  for(jp[1] = 0; jp[1] < 2; jp[1]++){
147  if(jp[0] == jp[1] && Nhem[!jp[0]] == 0) continue;
148  TLorentzVector hem_probes[2];
149  for(int i = 0; i < 2; i++) hem_probes[i] = hem[i];
150  for(int i = 0; i < 2; i++) hem_probes[jp[i]] += inputs[ip[i]];
151  double metric = hem_probes[0].P() + hem_probes[1].P();
152  if(metric > metric_max){
153  metric_max = metric;
154  for(int i = 0; i < 2; i++) ip_max[i] = ip[i];
155  for(int i = 0; i < 2; i++) jp_max[i] = jp[i];
156  }
157  }
158  }
159  }
160  }
161  if(metric_max < 0){
162  return false;
163  }
164 
165  // initialize output states
166  for(int i = 0; i < 2; i++) GetChildState(i).ClearElements();
167  for(int i = 0; i < 2; i++){
168  GetChildState(jp_max[i]).AddElement(GetInputState(ip_max[i]));
169  }
170  TVector3 nRef = inputs[ip_max[0]].Vect().Cross(inputs[ip_max[1]].Vect());
171  for(int i = 0; i < Ninput; i++){
172  if((i == ip_max[0]) || (i == ip_max[1])) continue;
173  int ihem = int(inputs[i].Vect().Dot(nRef) > 0.);
174  GetChildState(ihem).AddElement(GetInputState(i));
175  }
176  if(GetChildState(1).GetFourVector().M() > GetChildState(1).GetFourVector().M()){
177  std::vector<VisibleStateList> flip;
178  for(int i = 0; i < 2; i++) flip.push_back(GetChildState(i).GetElements());
179  for(int i = 0; i < 2; i++) GetChildState(i).ClearElements();
180  for(int i = 0; i < 2; i++) GetChildState(i).AddElements(flip[(i+1)%2]);
181  }
182 
183  ExecuteDependancyJigsaws();
184 
185  return SetSpirit(true);
186  }
187 
188  bool MinMassesCombJigsaw::EvaluateMetric(double& metric) const {
189  TLorentzVector P1 = GetDependancyStates(0).GetFourVector();
190  TLorentzVector P2 = GetDependancyStates(1).GetFourVector();
191 
192  double P = GetP((P1+P2).M(), P1.M(), P2.M());
193  if(P <= 0)
194  metric = -1.;
195  else
196  metric = 1./P;
197 
198  return true;
199  }
200 
201  MinMassesCombJigsaw MinMassesCombJigsaw::m_Empty;
202 
203 }
204