Fork me on GitHub

source: svn/trunk/modules/Weighter.cc@ 1379

Last change on this file since 1379 was 1125, checked in by Pavel Demin, 11 years ago

add module Weighter

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 3.9 KB
Line 
1
2/** \class Weighter
3 *
4 * Apply a weight depending on PDG code.
5 *
6 * $Date: 2013-05-26 22:36:00 +0000 (Sun, 26 May 2013) $
7 * $Revision: 1125 $
8 *
9 *
10 * \author P. Demin - UCL, Louvain-la-Neuve
11 *
12 */
13
14#include "modules/Weighter.h"
15
16#include "classes/DelphesClasses.h"
17#include "classes/DelphesFactory.h"
18#include "classes/DelphesFormula.h"
19
20#include "ExRootAnalysis/ExRootResult.h"
21#include "ExRootAnalysis/ExRootFilter.h"
22#include "ExRootAnalysis/ExRootClassifier.h"
23
24#include "TMath.h"
25#include "TString.h"
26#include "TFormula.h"
27#include "TRandom3.h"
28#include "TObjArray.h"
29#include "TDatabasePDG.h"
30#include "TLorentzVector.h"
31
32#include <algorithm>
33#include <stdexcept>
34#include <iostream>
35#include <sstream>
36
37using namespace std;
38
39//------------------------------------------------------------------------------
40
41bool Weighter::TIndexStruct::operator< (const Weighter::TIndexStruct &value) const
42{
43 Int_t i;
44
45 for(i = 0; i < 4; ++i)
46 {
47 if(codes[i] != value.codes[i]) return codes[i] < value.codes[i];
48 }
49
50 return false;
51}
52
53//------------------------------------------------------------------------------
54
55Weighter::Weighter() :
56 fItInputArray(0)
57{
58}
59
60//------------------------------------------------------------------------------
61
62Weighter::~Weighter()
63{
64}
65
66//------------------------------------------------------------------------------
67
68void Weighter::Init()
69{
70 ExRootConfParam param, paramCodes;
71 Int_t i, j, size, sizeCodes;
72 Int_t code;
73 TIndexStruct index;
74 Double_t weight;
75
76 fWeightSet.clear();
77
78
79 // set default weight value
80 fWeightMap.clear();
81 memset(index.codes, 0, 4*sizeof(Int_t));
82 fWeightMap[index] = 1.0;
83
84 // read weights
85 param = GetParam("Weight");
86 size = param.GetSize();
87 for(i = 0; i < size/2; ++i)
88 {
89 paramCodes = param[i*2];
90 sizeCodes = paramCodes.GetSize();
91 weight = param[i*2 + 1].GetDouble();
92
93 if(sizeCodes < 1 || sizeCodes > 4)
94 {
95 throw runtime_error("only 1, 2, 3 or 4 PDG codes can be specified per weight");
96 }
97
98 memset(index.codes, 0, 4*sizeof(Int_t));
99
100 for(j = 0; j < sizeCodes; ++j)
101 {
102 code = paramCodes[j].GetInt();
103 index.codes[j] = code;
104 fWeightSet.insert(code);
105 }
106
107 sort(index.codes, index.codes + 4);
108
109 fWeightMap[index] = weight;
110 }
111
112 // import input array(s)
113
114 fInputArray = ImportArray(GetString("InputArray", "Delphes/allParticles"));
115 fItInputArray = fInputArray->MakeIterator();
116
117 // create output array(s)
118
119 fOutputArray = ExportArray(GetString("OutputArray", "weight"));
120}
121
122//------------------------------------------------------------------------------
123
124void Weighter::Finish()
125{
126 if(fItInputArray) delete fItInputArray;
127}
128
129//------------------------------------------------------------------------------
130
131void Weighter::Process()
132{
133 Candidate *candidate;
134 Int_t i;
135 TIndexStruct index;
136 Double_t weight;
137 set<Int_t>::iterator itCodeSet;
138 map<TIndexStruct, Double_t>::iterator itWeightMap;
139
140 DelphesFactory *factory = GetFactory();
141
142 // loop over all particles
143 fCodeSet.clear();
144 fItInputArray->Reset();
145 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
146 {
147 if(candidate->Status != 3) continue;
148
149 if(fWeightSet.find(candidate->PID) == fWeightSet.end()) continue;
150
151 fCodeSet.insert(candidate->PID);
152 }
153
154 // find default weight value
155 memset(index.codes, 0, 4*sizeof(Int_t));
156 itWeightMap = fWeightMap.find(index);
157 weight = itWeightMap->second;
158
159 if(fCodeSet.size() <= 4)
160 {
161 i = 0;
162 for(itCodeSet = fCodeSet.begin(); itCodeSet != fCodeSet.end(); ++itCodeSet)
163 {
164 index.codes[i] = *itCodeSet;
165 ++i;
166 }
167
168 sort(index.codes, index.codes + 4);
169
170 itWeightMap = fWeightMap.find(index);
171 if(itWeightMap != fWeightMap.end())
172 {
173 weight = itWeightMap->second;
174 }
175 }
176
177 candidate = factory->NewCandidate();
178 candidate->Momentum.SetPtEtaPhiE(weight, 0.0, 0.0, weight);
179 fOutputArray->Add(candidate);
180}
181
182//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.