CbmRoot
PairAnalysisObjectCuts.cxx
Go to the documentation of this file.
1 //
3 //
4 // Authors:
5 // Julian Book <Julian.Book@cern.ch>
6 /*
7 
8  Advanced cut class.
9 
10  Add cuts that depend on some variable or formula of them. The minimum and
11  maximum cut value can be expressed by a TForumla, TGraph or a THnBase.
12 
13 
14 */
15 // //
17 
18 
19 #include <TAxis.h>
20 #include <TFormula.h>
21 #include <TGraph.h>
22 #include <THnBase.h>
23 #include <TSpline.h>
24 
25 #include "PairAnalysisHelper.h"
26 #include "PairAnalysisObjectCuts.h"
27 
28 ClassImp(PairAnalysisObjectCuts)
29 
30 
31  PairAnalysisObjectCuts::PairAnalysisObjectCuts()
32  : PairAnalysisObjectCuts("objcuts", "objcuts") {
33  //
34  // Default costructor
35  //
36 }
37 
38 //________________________________________________________________________
39 PairAnalysisObjectCuts::PairAnalysisObjectCuts(const char* name,
40  const char* title)
41  : AnalysisCuts(name, title)
42  , fUsedVars(new TBits(PairAnalysisVarManager::kNMaxValuesMC)) {
43  //
44  // Named contructor
45  //
46  for (Int_t i = 0; i < PairAnalysisObjectCuts::kNMaxCuts; ++i) {
47  fActiveCuts[i] = 0;
48  fCutExclude[i] = kFALSE;
49  fCutMin[i] = 0x0;
50  fCutMax[i] = 0x0;
51  fVarFormula[i] = 0x0;
52  }
53  PairAnalysisVarManager::InitFormulas();
54 }
55 
56 //________________________________________________________________________
57 PairAnalysisObjectCuts::~PairAnalysisObjectCuts() {
58  //
59  // Destructor
60  //
61  if (fUsedVars) delete fUsedVars;
62  for (Int_t i = 0; i < PairAnalysisObjectCuts::kNMaxCuts; ++i) {
63  fActiveCuts[i] = 0;
64  if (fCutMin[i]) delete fCutMin[i];
65  if (fCutMax[i]) delete fCutMax[i];
66  if (fVarFormula[i]) delete fVarFormula[i];
67  }
68 }
69 
70 //________________________________________________________________________
71 Bool_t PairAnalysisObjectCuts::IsSelected(Double_t* const values) {
72  //
73  // Make cut decision
74  //
75 
76  //reset
77  fSelectedCutsMask = 0;
78  SetSelected(kFALSE);
79 
80  for (Int_t iCut = 0; iCut < fNActiveCuts; ++iCut) {
81  Int_t cut = fActiveCuts[iCut];
82  SETBIT(fSelectedCutsMask, iCut);
83 
84  // variable you want to cut on
85  Double_t compValue = 0.;
86  if (fVarFormula[iCut])
87  compValue = PairAnalysisHelper::EvalFormula(fVarFormula[iCut], values);
88  else
89  compValue = values[cut];
90 
91  // cut limits
92  Double_t cutMin = -9.e30;
93  Double_t cutMax = +9.e30;
94 
95 
97  if ((fCutMin[iCut] && fCutMin[iCut]->IsA() == THnBase::Class())
98  || (fCutMax[iCut] && fCutMax[iCut]->IsA() == THnBase::Class())) {
99 
100  THnBase* hn = static_cast<THnBase*>(fCutMin[iCut]);
101  if (fCutMin[iCut]) {
102  Double_t* vals = new Double_t[hn->GetNdimensions()]; //={-1};
103  // get array of values for the corresponding dimensions using axis names
104  for (Int_t idim = 0; idim < hn->GetNdimensions(); idim++) {
105  vals[idim] = values[PairAnalysisVarManager::GetValueType(
106  hn->GetAxis(idim)
107  ->GetName())];
108  }
109  // find bin for values (w/o creating it in case it is not filled)
110  Long_t bin = hn->GetBin(vals, kFALSE);
111  if (bin > 0) cutMin = hn->GetBinContent(bin);
112  delete[] vals;
113  }
114 
115  THnBase* hx = static_cast<THnBase*>(fCutMax[iCut]);
116  if (fCutMax[iCut]) {
117  Double_t* vals = new Double_t[hx->GetNdimensions()]; //={-1};
118  // get array of values for the corresponding dimensions using axis names
119  for (Int_t idim = 0; idim < hx->GetNdimensions(); idim++) {
120  vals[idim] = values[PairAnalysisVarManager::GetValueType(
121  hx->GetAxis(idim)
122  ->GetName())];
123  }
124  // find bin for values (w/o creating it in case it is not filled)
125  Long_t bin = hx->GetBin(vals, kFALSE);
126  if (bin > 0) cutMax = hx->GetBinContent(bin);
127  delete[] vals;
128  }
129 
130  } else if ((fCutMin[iCut] && fCutMin[iCut]->IsA() == TFormula::Class())
131  || (fCutMax[iCut]
132  && fCutMax[iCut]->IsA() == TFormula::Class())) {
134  TFormula* formN = static_cast<TFormula*>(fCutMin[iCut]);
135  if (formN) cutMin = PairAnalysisHelper::EvalFormula(formN, values);
136 
137  TFormula* formM = static_cast<TFormula*>(fCutMax[iCut]);
138  if (formM) cutMax = PairAnalysisHelper::EvalFormula(formM, values);
139 
140  } else if ((fCutMin[iCut] && fCutMin[iCut]->IsA() == TGraph::Class())
141  || (fCutMax[iCut] && fCutMax[iCut]->IsA() == TGraph::Class())) {
145  TGraph* graphN = static_cast<TGraph*>(fCutMin[iCut]);
146  TGraph* graphM = static_cast<TGraph*>(fCutMax[iCut]);
147 
148  // get x-value from formula
149  TFormula* formX = 0;
150  Double_t xval = 0.;
151  if (graphN) {
152  formX = static_cast<TFormula*>(graphN->GetListOfFunctions()->At(0));
153  if (formX) xval = PairAnalysisHelper::EvalFormula(formX, values);
154  }
155  if (!formX && graphM) {
156  formX = static_cast<TFormula*>(graphM->GetListOfFunctions()->At(0));
157  if (formX) xval = PairAnalysisHelper::EvalFormula(formX, values);
158  }
159 
161  // if(graphN) cutMin = graphN->Eval(xval);
162  // if(graphM) cutMax = graphM->Eval(xval);
163 
165  Int_t idx = TMath::BinarySearch(graphN->GetN(), graphN->GetX(), xval);
166  cutMin = graphN->GetY()[idx];
167  idx = TMath::BinarySearch(graphM->GetN(), graphM->GetX(), xval);
168  cutMax = graphM->GetY()[idx];
169 
170  } else if ((fCutMin[iCut] && fCutMin[iCut]->IsA() == TSpline3::Class())
171  || (fCutMax[iCut]
172  && fCutMax[iCut]->IsA() == TSpline3::Class())) {
175  // TSpline3 *splineN = static_cast<TSpline3*>(fCutMin[iCut]);
176  // if(splineN) cutMin = splineN->Eval(xval);
177  // TSpline3 *splineM = static_cast<TSpline3*>(fCutMax[iCut]);
178  // if(splineM) cutMax = splineM->Eval(xval);
179  } else {
180  Error("IsSelected:",
181  "Cut object not supported (this message should never appear)");
182  return kTRUE;
183  }
184 
185  // protection against NaN (e.g. outside formula ranges, exclude them)
186  if (TMath::IsNaN(cutMin)) cutMin = compValue + 1.;
187  if (TMath::IsNaN(cutMax)) cutMax = compValue - 1.;
188 
189  // apply cut
190  if (((compValue < cutMin) || (compValue > cutMax)) ^ fCutExclude[iCut])
191  CLRBIT(fSelectedCutsMask, iCut);
192 
193  // cut type and decision
194  if (fCutType == kAll && !TESTBIT(fSelectedCutsMask, iCut))
195  return kFALSE; // option to (minor) speed improvement
196  }
197 
198  Bool_t isSelected = (fSelectedCutsMask == fActiveCutsMask);
199  if (fCutType == kAny) isSelected = (fSelectedCutsMask > 0);
200  SetSelected(isSelected);
201  return isSelected;
202 }
203 
204 //________________________________________________________________________
205 Bool_t PairAnalysisObjectCuts::IsSelected(TObject* track) {
206  //
207  // Make cut decision
208  //
209 
210  //reset
211  fSelectedCutsMask = 0;
212  SetSelected(kFALSE);
213  if (!track) return kFALSE;
214 
215  //Fill values
216  Double_t* values = PairAnalysisVarManager::GetData();
217  PairAnalysisVarManager::SetFillMap(fUsedVars);
218  PairAnalysisVarManager::Fill(track, values);
219 
220  return (IsSelected(values));
221 }
222 
223 //________________________________________________________________________
224 void PairAnalysisObjectCuts::AddCut(PairAnalysisVarManager::ValueTypes type,
225  const char* formulaMin,
226  const char* formulaMax,
227  Bool_t excludeRange) {
228  //
229  // Set cut range and activate it
230  //
231  fCutMin[fNActiveCuts] =
232  PairAnalysisHelper::GetFormula(formulaMin, formulaMin);
233  fCutMax[fNActiveCuts] =
234  PairAnalysisHelper::GetFormula(formulaMax, formulaMax);
235  fCutExclude[fNActiveCuts] = excludeRange;
236  SETBIT(fActiveCutsMask, fNActiveCuts);
237  fActiveCuts[fNActiveCuts] = (UShort_t) type;
238  fUsedVars->SetBitNumber(type, kTRUE);
239  ++fNActiveCuts;
240 }
241 
242 //________________________________________________________________________
243 void PairAnalysisObjectCuts::AddCut(const char* formula,
244  const char* formulaMin,
245  const char* formulaMax,
246  Bool_t excludeRange) {
247  //
248  // Set cut range and activate it
249  //
250  fCutMin[fNActiveCuts] =
251  PairAnalysisHelper::GetFormula(formulaMin, formulaMin);
252  fCutMax[fNActiveCuts] =
253  PairAnalysisHelper::GetFormula(formulaMax, formulaMax);
254  fCutExclude[fNActiveCuts] = excludeRange;
255  SETBIT(fActiveCutsMask, fNActiveCuts);
256  fActiveCuts[fNActiveCuts] = (UShort_t) PairAnalysisVarManager::kNMaxValuesMC;
257  fVarFormula[fNActiveCuts] = PairAnalysisHelper::GetFormula(formula, formula);
258  ++fNActiveCuts;
259 }
260 
261 //________________________________________________________________________
262 void PairAnalysisObjectCuts::AddCut(PairAnalysisVarManager::ValueTypes type,
263  TGraph* const graphMin,
264  TGraph* const graphMax,
265  Bool_t excludeRange) {
266  //
267  // Set cut range and activate it
268  //
269  // if(graphMin) {
270  // TSpline3 *splineMin = new TSpline3("splineMin",graphMin); //spline w/o begin and end point conditions
271  // fCutMin[fNActiveCuts]=splineMin;
272  // }
273  // if(graphMax) {
274  // TSpline3 *splineMax = new TSpline3("splineMax",graphMax); //spline w/o begin and end point conditions
275  // fCutMax[fNActiveCuts]=splineMax;
276  // }
277  fCutMin[fNActiveCuts] = graphMin;
278  fCutMax[fNActiveCuts] = graphMax;
279  fCutExclude[fNActiveCuts] = excludeRange;
280  SETBIT(fActiveCutsMask, fNActiveCuts);
281  fActiveCuts[fNActiveCuts] = (UShort_t) type;
282  fUsedVars->SetBitNumber(type, kTRUE);
283  ++fNActiveCuts;
284 }
285 
286 //________________________________________________________________________
287 void PairAnalysisObjectCuts::AddCut(const char* formula,
288  TGraph* const graphMin,
289  TGraph* const graphMax,
290  Bool_t excludeRange) {
291  //
292  // Set cut range and activate it
293  //
294  // if(graphMin) {
295  // TSpline3 *splineMin = new TSpline3("splineMin",graphMin); //spline w/o begin and end point conditions
296  // fCutMin[fNActiveCuts]=splineMin;
297  // }
298  // if(graphMax) {
299  // TSpline3 *splineMax = new TSpline3("splineMax",graphMax); //spline w/o begin and end point conditions
300  // fCutMax[fNActiveCuts]=splineMax;
301  // }
302  fCutMin[fNActiveCuts] = graphMin;
303  fCutMax[fNActiveCuts] = graphMax;
304  fCutExclude[fNActiveCuts] = excludeRange;
305  SETBIT(fActiveCutsMask, fNActiveCuts);
306  fActiveCuts[fNActiveCuts] = (UShort_t) PairAnalysisVarManager::kNMaxValuesMC;
307  fVarFormula[fNActiveCuts] = PairAnalysisHelper::GetFormula(formula, formula);
308  ++fNActiveCuts;
309 }
310 
311 //________________________________________________________________________
312 void PairAnalysisObjectCuts::AddCut(PairAnalysisVarManager::ValueTypes type,
313  THnBase* const histMin,
314  THnBase* const histMax,
315  Bool_t excludeRange) {
316  //
317  // Set cut range and activate it
318  //
319  fCutMin[fNActiveCuts] = histMin;
320  fCutMax[fNActiveCuts] = histMax;
321  fCutExclude[fNActiveCuts] = excludeRange;
322  SETBIT(fActiveCutsMask, fNActiveCuts);
323  fActiveCuts[fNActiveCuts] = (UShort_t) type;
324  fUsedVars->SetBitNumber(type, kTRUE);
325  ++fNActiveCuts;
326 }
327 
328 //________________________________________________________________________
329 void PairAnalysisObjectCuts::AddCut(const char* formula,
330  THnBase* const histMin,
331  THnBase* const histMax,
332  Bool_t excludeRange) {
333  //
334  // Set cut range and activate it
335  //
336  fCutMin[fNActiveCuts] = histMin;
337  fCutMax[fNActiveCuts] = histMax;
338  fCutExclude[fNActiveCuts] = excludeRange;
339  SETBIT(fActiveCutsMask, fNActiveCuts);
340  fActiveCuts[fNActiveCuts] = (UShort_t) PairAnalysisVarManager::kNMaxValuesMC;
341  fVarFormula[fNActiveCuts] = PairAnalysisHelper::GetFormula(formula, formula);
342  ++fNActiveCuts;
343 }
344 
345 //________________________________________________________________________
346 void PairAnalysisObjectCuts::Print(const Option_t* /*option*/) const {
347  //
348  // Print cuts and the range
349  //
350  printf("cut ranges for '%s'\n", GetTitle());
351  if (fCutType == kAll) {
352  printf("All Cuts have to be fulfilled\n");
353  } else {
354  printf("Any Cut has to be fulfilled\n");
355  }
356 
358  for (Int_t iCut = 0; iCut < fNActiveCuts; ++iCut) {
359 
360  // variable
361  Int_t cut = (Int_t) fActiveCuts[iCut];
362  TString tit = PairAnalysisVarManager::GetValueName((Int_t) cut);
363 
364  Bool_t fvar = fVarFormula[iCut];
365  if (fvar) {
366  TFormula* form = static_cast<TFormula*>(fVarFormula[iCut]);
367  tit = form->GetExpFormula();
368  // replace parameter variables with names labels
369  for (Int_t j = 0; j < form->GetNpar(); j++)
370  tit.ReplaceAll(Form("[%d]", j), form->GetParName(j));
371  }
372 
373  // cut logic
374  Bool_t inverse = fCutExclude[iCut];
375 
376  // cut limits
377  Bool_t bCutGraph =
378  (fCutMin[iCut] && fCutMin[iCut]->IsA() == TGraph::Class());
379  Bool_t bCutForm =
380  (fCutMin[iCut] && fCutMin[iCut]->IsA() == TFormula::Class());
381  Bool_t bCutHn = (fCutMin[iCut] && fCutMin[iCut]->IsA() == THnBase::Class());
382  // Bool_t bCutSpline = (fCutMin[iCut] && fCutMin[iCut]->IsA() == TSpline::Class());
383 
384  TString dep = "";
385  // HnBase
386  if (bCutHn) {
387  THnBase* obj = static_cast<THnBase*>(fCutMin[iCut]);
388  for (Int_t idim = 0; idim < obj->GetNdimensions(); idim++)
389  dep += Form("%s%s", (idim ? "," : ""), obj->GetAxis(idim)->GetName());
390  dep.Prepend("histogram(");
391  dep.Append(")");
392  }
393  // Graph
394  if (bCutGraph) {
395  TGraph* obj = static_cast<TGraph*>(fCutMin[iCut]);
396  TFormula* form = static_cast<TFormula*>(obj->GetListOfFunctions()->At(0));
397  dep = form->GetExpFormula();
398  // replace parameter variables with names labels
399  for (Int_t j = 0; j < form->GetNpar(); j++)
400  dep.ReplaceAll(Form("[%d]", j), form->GetParName(j));
401  dep.Prepend("graph(");
402  dep.Append(")");
403  }
404  // formula
405  if (bCutForm) {
406  TFormula* obj = static_cast<TFormula*>(fCutMin[iCut]);
407  dep = obj->GetExpFormula();
408  // replace parameter variables with names labels
409  for (Int_t j = 0; j < obj->GetNpar(); j++)
410  dep.ReplaceAll(Form("[%d]", j), obj->GetParName(j));
411  dep.Prepend("formula(");
412  dep.Append(")");
413  }
414 
415  // stdout
416  if (!inverse)
417  printf(
418  "Cut %02d: %s < %s < %s\n", iCut, dep.Data(), tit.Data(), dep.Data());
419  else
420  printf("Cut %02d: !(%s < %s < %s)\n",
421  iCut,
422  dep.Data(),
423  tit.Data(),
424  dep.Data());
425 
426  } //loop over cuts
427 }
PairAnalysisObjectCuts::fCutMax
TObject * fCutMax[PairAnalysisObjectCuts::kNMaxCuts]
Definition: PairAnalysisObjectCuts.h:97
PairAnalysisObjectCuts::fSelectedCutsMask
UInt_t fSelectedCutsMask
Definition: PairAnalysisObjectCuts.h:89
PairAnalysisVarManager::kNMaxValuesMC
@ kNMaxValuesMC
Definition: PairAnalysisVarManager.h:343
PairAnalysisVarManager::SetFillMap
static void SetFillMap(TBits *map)
Definition: PairAnalysisVarManager.h:368
PairAnalysisHelper::EvalFormula
Double_t EvalFormula(TFormula *form, const Double_t *values)
Definition: PairAnalysisHelper.cxx:278
PairAnalysisObjectCuts::IsSelected
virtual Bool_t IsSelected(Double_t *const values)
Definition: PairAnalysisObjectCuts.cxx:71
PairAnalysisObjectCuts::fVarFormula
TFormula * fVarFormula[PairAnalysisObjectCuts::kNMaxCuts]
Definition: PairAnalysisObjectCuts.h:99
PairAnalysisHelper.h
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
PairAnalysisVarManager
Definition: PairAnalysisVarManager.h:68
PairAnalysisVarManager::Fill
static void Fill(const TObject *particle, Double_t *const values)
Definition: PairAnalysisVarManager.h:474
PairAnalysisVarManager::GetValueName
static const char * GetValueName(Int_t i)
Definition: PairAnalysisVarManager.h:377
PairAnalysisVarManager::GetValueType
static UInt_t GetValueType(const char *valname)
Definition: PairAnalysisVarManager.cxx:355
PairAnalysisObjectCuts::kAll
@ kAll
Definition: PairAnalysisObjectCuts.h:25
PairAnalysisObjectCuts::fCutExclude
Bool_t fCutExclude[PairAnalysisObjectCuts::kNMaxCuts]
Definition: PairAnalysisObjectCuts.h:93
PairAnalysisVarManager::InitFormulas
static void InitFormulas()
Definition: PairAnalysisVarManager.h:1692
PairAnalysisObjectCuts::~PairAnalysisObjectCuts
virtual ~PairAnalysisObjectCuts()
Definition: PairAnalysisObjectCuts.cxx:57
PairAnalysisHelper::GetFormula
TFormula * GetFormula(const char *name, const char *formula)
Definition: PairAnalysisHelper.cxx:330
PairAnalysisObjectCuts::fActiveCutsMask
UInt_t fActiveCutsMask
Definition: PairAnalysisObjectCuts.h:87
AnalysisCuts::SetSelected
virtual void SetSelected(Bool_t dec)
Definition: AnalysisCuts.h:25
ClassImp
ClassImp(PairAnalysisObjectCuts) PairAnalysisObjectCuts
Definition: PairAnalysisObjectCuts.cxx:28
PairAnalysisObjectCuts::kAny
@ kAny
Definition: PairAnalysisObjectCuts.h:25
PairAnalysisObjectCuts::fCutMin
TObject * fCutMin[PairAnalysisObjectCuts::kNMaxCuts]
Definition: PairAnalysisObjectCuts.h:95
PairAnalysisObjectCuts
Definition: PairAnalysisObjectCuts.h:22
PairAnalysisObjectCuts::fUsedVars
TBits * fUsedVars
Definition: PairAnalysisObjectCuts.h:83
PairAnalysisObjectCuts::fActiveCuts
UShort_t fActiveCuts[PairAnalysisObjectCuts::kNMaxCuts]
Definition: PairAnalysisObjectCuts.h:85
PairAnalysisVarManager::GetData
static Double_t * GetData()
Definition: PairAnalysisVarManager.h:386
PairAnalysisObjectCuts.h
PairAnalysisObjectCuts::PairAnalysisObjectCuts
PairAnalysisObjectCuts()
PairAnalysisObjectCuts::fCutType
CutType fCutType
Definition: PairAnalysisObjectCuts.h:91
PairAnalysisVarManager::ValueTypes
ValueTypes
Definition: PairAnalysisVarManager.h:71
PairAnalysisObjectCuts::Print
virtual void Print(const Option_t *option="") const
Definition: PairAnalysisObjectCuts.cxx:346
PairAnalysisObjectCuts::AddCut
void AddCut(PairAnalysisVarManager::ValueTypes type, const char *formulaMin, const char *formulaMax, Bool_t excludeRange=kFALSE)
Definition: PairAnalysisObjectCuts.cxx:224
AnalysisCuts
Definition: AnalysisCuts.h:12
PairAnalysisObjectCuts::fNActiveCuts
UShort_t fNActiveCuts
Definition: PairAnalysisObjectCuts.h:86
PairAnalysisObjectCuts::kNMaxCuts
@ kNMaxCuts
Definition: PairAnalysisObjectCuts.h:26