ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitAna/PhysicsMod/src/FwdJetEvtSelMod.cc
Revision: 1.2
Committed: Thu Oct 9 18:10:23 2008 UTC (16 years, 7 months ago) by ceballos
Content type: text/plain
Branch: MAIN
Changes since 1.1: +6 -3 lines
Log Message: 
Fixing memory leaks here and there

File Contents

# User Rev Content
1 ceballos 1.2 // $Id: FwdJetEvtSelMod.cc,v 1.1 2008/10/06 16:59:48 ceballos Exp $
2 ceballos 1.1
3     #include "MitAna/PhysicsMod/interface/FwdJetEvtSelMod.h"
4     #include <TH1D.h>
5     #include <TH2D.h>
6     #include "MitAna/DataTree/interface/Names.h"
7     #include "MitAna/DataCont/interface/ObjArray.h"
8     #include "MitCommon/MathTools/interface/MathUtils.h"
9     #include "MitAna/Utils/interface/DiTauSystem.h"
10    
11     using namespace mithep;
12     ClassImp(mithep::FwdJetEvtSelMod)
13    
14     //--------------------------------------------------------------------------------------------------
15     FwdJetEvtSelMod::FwdJetEvtSelMod(const char *name, const char *title) :
16     BaseMod(name,title),
17     fPrintDebug(false),
18     fFillHist(false),
19     fMetName(Names::gkCaloMetBrn),
20     fCleanJetsName(Names::gkCleanJetsName),
21     fCleanFwdJetsName(Names::gkCleanFwdJetsName),
22     fCleanNoFwdJetsName(Names::gkCleanNoFwdJetsName),
23     fMCqqHsName(Names::gkMCqqHsName),
24     fMet(0),
25     fUseANN(false),
26     fUseHighestPtJets(true),
27     fJetPtMax(40),
28     fJetPtMin(40),
29     fDeltaEtaMin(4),
30     fDiJetMassMin(900),
31     fANNOutputMin(0.5),
32     fNEventsProcessed(0)
33     {
34     // Constructor.
35     }
36    
37     //-----------------------------------------------------------------------------
38     void FwdJetEvtSelMod::Begin()
39     {
40     // Run startup code on the client machine. For this module, we dont do
41     // anything here.
42     }
43    
44     //-----------------------------------------------------------------------------
45     void FwdJetEvtSelMod::Process()
46     {
47     // Process entries of the tree. For this module, we just load the branches and
48     fNEventsProcessed++;
49    
50     if (fNEventsProcessed % 1000 == 0 || fPrintDebug) {
51     time_t systime;
52     systime = time(NULL);
53     cerr << endl << "FwdJetEvtSelMod : Process Event " << fNEventsProcessed << " Time: " << ctime(&systime) << endl;
54     }
55    
56     //Get Generator Level information for matching
57     ObjArray<MCParticle> *GenqqHs = dynamic_cast<ObjArray<MCParticle>* >
58     (FindObjThisEvt(fMCqqHsName.Data()));
59    
60     //Obtain all the good objects from the event cleaning module
61     ObjArray<Electron> *CleanElectrons = dynamic_cast<ObjArray<Electron>* >
62     (FindObjThisEvt(Names::gkCleanElectronsName));
63     ObjArray<Muon> *CleanMuons = dynamic_cast<ObjArray<Muon>* >
64     (FindObjThisEvt(Names::gkCleanMuonsName));
65     ObjArray<Jet> *CleanJets = dynamic_cast<ObjArray<Jet>* >
66     (FindObjThisEvt(fCleanJetsName.Data()));
67    
68 ceballos 1.2 ObjArray<Jet> *CleanFwdJets = new ObjArray<Jet>; CleanFwdJets->SetOwner(true);
69     ObjArray<Jet> *CleanNoFwdJets = new ObjArray<Jet>; CleanNoFwdJets->SetOwner(true);
70 ceballos 1.1
71     // Sort and count the number of central Jets for vetoing
72     vector<Jet*> sortedJets;
73     for(UInt_t i=0; i<CleanJets->GetEntries(); i++){
74     Jet* jet_f = new Jet(CleanJets->At(i)->Px()*CleanJets->At(i)->L2RelativeCorrectionScale()*CleanJets->At(i)->L3AbsoluteCorrectionScale(),
75     CleanJets->At(i)->Py()*CleanJets->At(i)->L2RelativeCorrectionScale()*CleanJets->At(i)->L3AbsoluteCorrectionScale(),
76     CleanJets->At(i)->Pz()*CleanJets->At(i)->L2RelativeCorrectionScale()*CleanJets->At(i)->L3AbsoluteCorrectionScale(),
77     CleanJets->At(i)->E() *CleanJets->At(i)->L2RelativeCorrectionScale()*CleanJets->At(i)->L3AbsoluteCorrectionScale());
78     sortedJets.push_back(jet_f);
79     }
80     for(UInt_t i=0; i<sortedJets.size(); i++){
81     for(UInt_t j=i+1; j<sortedJets.size(); j++){
82     if(sortedJets[i]->Pt() < sortedJets[j]->Pt()) {
83     //swap i and j
84     Jet* tempjet = sortedJets[i];
85     sortedJets[i] = sortedJets[j];
86     sortedJets[j] = tempjet;
87     }
88     }
89     }
90     if(sortedJets.size() >= 2 && sortedJets[0]->Pt() > 20 &&
91     sortedJets[1]->Pt() > 20){
92     // Use either the 2 or the 4 highest Pt jets
93     UInt_t NMaxLoop = 4;
94     if(fUseHighestPtJets == true) NMaxLoop = 2;
95     int jmax = -1;
96     int jmin = -1;
97     double nnOutputPair = 0.0;
98     for(UInt_t i=0; i<TMath::Min(sortedJets.size()-1,NMaxLoop-1); i++){
99     for(UInt_t j=i+1; j<TMath::Min(sortedJets.size(),NMaxLoop); j++){
100     // Stop if a good jet pair is already found
101     if(jmax != -1) continue;
102     CompositeParticle *dijet = new CompositeParticle();
103     dijet->AddDaughter(sortedJets[i]);
104     dijet->AddDaughter(sortedJets[j]);
105     double deltaPhi = MathUtils::DeltaPhi(sortedJets[i]->Phi(),
106     sortedJets[j]->Phi());
107     double deltaEta = fabs(sortedJets[i]->Eta()-sortedJets[j]->Eta());
108     bool isFwdJet[2] = {false, false};
109     isFwdJet[0] = sortedJets[i]->Pt() > fJetPtMax &&
110     sortedJets[j]->Pt() > fJetPtMin &&
111     sortedJets[i]->Eta()*sortedJets[j]->Eta() < 0 &&
112     deltaEta > fDeltaEtaMin && dijet->Mass() > fDiJetMassMin;
113     double varqqHNN[6] = {sortedJets[i]->Pt(),sortedJets[j]->Pt(),
114     deltaEta,deltaPhi,dijet->Mass(),
115     sortedJets[i]->Eta()*sortedJets[j]->Eta()};
116     if(deltaEta > 2.0 && dijet->Mass() > 300)
117     nnOutputPair = Testmlp_qqH(varqqHNN);
118     isFwdJet[1] = nnOutputPair > fANNOutputMin;
119     if (fUseANN == false and isFwdJet[0] == true){
120     CleanFwdJets->Add(sortedJets[i]);
121     CleanFwdJets->Add(sortedJets[j]);
122     jmax = i; jmin = j;
123     }
124     else if(fUseANN == true and isFwdJet[1] == true){
125     CleanFwdJets->Add(sortedJets[i]);
126     CleanFwdJets->Add(sortedJets[j]);
127     jmax = i; jmin = j;
128     }
129 ceballos 1.2 delete dijet;
130 ceballos 1.1 } // j...
131     } // i...
132     // Looking for central jets if a good forward tagging pair is found
133     if(jmax >= 0){
134     for(UInt_t i=0; i<sortedJets.size(); i++){
135     if((sortedJets[i]->Eta() > sortedJets[jmax]->Eta() &&
136     sortedJets[i]->Eta() < sortedJets[jmin]->Eta()) ||
137     (sortedJets[i]->Eta() > sortedJets[jmin]->Eta() &&
138     sortedJets[i]->Eta() < sortedJets[jmax]->Eta())){
139     CleanNoFwdJets->Add(sortedJets[i]);
140     }
141     }
142     }
143     if(fFillHist == true){
144     if(fUseHighestPtJets == true){
145     CompositeParticle *dijet = new CompositeParticle();
146     dijet->AddDaughter(sortedJets[0]);
147     dijet->AddDaughter(sortedJets[1]);
148     double deltaPhi = MathUtils::DeltaPhi(sortedJets[0]->Phi(),
149     sortedJets[1]->Phi());
150     double deltaEta = fabs(sortedJets[0]->Eta()-sortedJets[1]->Eta());
151     bool isFwdJet[2] = {false, false};
152     isFwdJet[0] = sortedJets[0]->Pt() > fJetPtMax &&
153     sortedJets[1]->Pt() > fJetPtMin &&
154     sortedJets[0]->Eta()*sortedJets[1]->Eta() < 0 &&
155     deltaEta > fDeltaEtaMin && dijet->Mass() > fDiJetMassMin;
156     double varqqHNN[6] = {sortedJets[0]->Pt(),sortedJets[1]->Pt(),
157     deltaEta,deltaPhi,dijet->Mass(),
158     sortedJets[0]->Eta()*sortedJets[1]->Eta()};
159     double nnOutput = 0.0;
160     if(deltaEta > 2.0 && dijet->Mass() > 300)
161     nnOutput = Testmlp_qqH(varqqHNN);
162     isFwdJet[1] = nnOutput > fANNOutputMin;
163     int VBYtype = 2;
164     if(GenqqHs->GetEntries() == 2){
165     double deltaR00 = MathUtils::DeltaR(sortedJets[0]->Phi(),sortedJets[0]->Eta(),
166     GenqqHs->At(0)->Phi(),GenqqHs->At(0)->Eta());
167     double deltaR01 = MathUtils::DeltaR(sortedJets[0]->Phi(),sortedJets[0]->Eta(),
168     GenqqHs->At(1)->Phi(),GenqqHs->At(1)->Eta());
169     double deltaR10 = MathUtils::DeltaR(sortedJets[1]->Phi(),sortedJets[1]->Eta(),
170     GenqqHs->At(0)->Phi(),GenqqHs->At(0)->Eta());
171     double deltaR11 = MathUtils::DeltaR(sortedJets[1]->Phi(),sortedJets[1]->Eta(),
172     GenqqHs->At(1)->Phi(),GenqqHs->At(1)->Eta());
173     if((deltaR00 < 0.5 && deltaR11 < 0.5) ||
174     (deltaR01 < 0.5 && deltaR10 < 0.5)) VBYtype = 0;
175     else VBYtype = 1;
176     }
177     hDFwdJetSel[0+VBYtype*100]->Fill(deltaPhi*180./TMath::Pi());
178     hDFwdJetSel[1+VBYtype*100]->Fill(TMath::Min(deltaEta,9.999));
179     hDFwdJetSel[2+VBYtype*100]->Fill(TMath::Min(sortedJets[0]->Pt(),399.999));
180     hDFwdJetSel[3+VBYtype*100]->Fill(TMath::Min(sortedJets[1]->Pt(),399.999));
181     hDFwdJetSel[4+VBYtype*100]->Fill(TMath::Min(dijet->Mass(),3999.999));
182     hDFwdJetSel[5+VBYtype*100]->Fill(sortedJets[0]->Eta()* sortedJets[1]->Eta()/
183     fabs(sortedJets[0]->Eta()* sortedJets[1]->Eta()));
184     hDFwdJetSel[6+VBYtype*100]->Fill(nnOutput);
185     hDFwdJetSel[7+VBYtype*100]->Fill((double)isFwdJet[0]);
186     hDFwdJetSel[8+VBYtype*100]->Fill((double)isFwdJet[1]);
187 ceballos 1.2 delete dijet;
188 ceballos 1.1 }
189    
190     // Event has a good forward tagging pair
191     if(jmax >= 0){
192     int VBYtype = 2;
193     if(GenqqHs->GetEntries() == 2){
194     double deltaR00 = MathUtils::DeltaR(sortedJets[jmax]->Phi(),sortedJets[jmax]->Eta(),
195     GenqqHs->At(0)->Phi(),GenqqHs->At(0)->Eta());
196     double deltaR01 = MathUtils::DeltaR(sortedJets[jmax]->Phi(),sortedJets[jmax]->Eta(),
197     GenqqHs->At(1)->Phi(),GenqqHs->At(1)->Eta());
198     double deltaR10 = MathUtils::DeltaR(sortedJets[jmin]->Phi(),sortedJets[jmin]->Eta(),
199     GenqqHs->At(0)->Phi(),GenqqHs->At(0)->Eta());
200     double deltaR11 = MathUtils::DeltaR(sortedJets[jmin]->Phi(),sortedJets[jmin]->Eta(),
201     GenqqHs->At(1)->Phi(),GenqqHs->At(1)->Eta());
202     if((deltaR00 < 0.5 && deltaR11 < 0.5) ||
203     (deltaR01 < 0.5 && deltaR10 < 0.5)) VBYtype = 0;
204     else VBYtype = 1;
205     }
206     hDFwdJetSel[9+VBYtype*100]->Fill((double)CleanNoFwdJets->GetEntries());
207     // Lepton selection
208     vector<ChargedParticle*> leptons;
209     vector<string> leptonType;
210    
211     // Make lepton vector from muons and electrons
212     for (UInt_t j=0; j<CleanMuons->GetEntries(); j++) {
213     leptons.push_back(CleanMuons->At(j));
214     leptonType.push_back("mu");
215     }
216     for (UInt_t j=0; j<CleanElectrons->GetEntries(); j++) {
217     leptons.push_back(CleanElectrons->At(j));
218     leptonType.push_back("e");
219     }
220     hDFwdJetSel[10+VBYtype*100]->Fill((double)leptons.size());
221    
222     if(leptons.size() >= 2){
223     // Sort the Leptons by Pt
224     for(UInt_t i=0; i<leptons.size(); i++){
225     for(UInt_t j=i+1; j<leptons.size(); j++){
226     if(leptons[i]->Pt() < leptons[j]->Pt()) {
227     //swap i and j
228     ChargedParticle* templepton = leptons[i];
229     leptons[i] = leptons[j];
230     leptons[j] = templepton;
231     string templeptonType = leptonType[i];
232     leptonType[i] = leptonType[j];
233     leptonType[j] = templeptonType;
234     }
235     }
236     }
237     CompositeParticle *dilepton = new CompositeParticle();
238     dilepton->AddDaughter(leptons[0]);
239     dilepton->AddDaughter(leptons[1]);
240    
241     int pairType = -1;
242     if (leptonType[0] == "mu" && leptonType[1] == "mu" )
243     pairType = 0;
244     else if(leptonType[0] == "e" && leptonType[1] == "e")
245     pairType = 1;
246     else if((leptonType[0] == "e" && leptonType[1] == "mu") ||
247     (leptonType[0] == "mu" && leptonType[1] == "e"))
248     pairType = 2;
249     else {
250     cout << "Hey, this is not possible, leptonTypes: "
251     << leptonType[0] << " - "
252     << leptonType[1] << endl;
253     }
254    
255     // Charge of the leptons should be opposite
256     if (dilepton->Charge() == 0 && dilepton->Mass() > 12 &&
257     leptons[0]->Pt() > 20 && leptons[0]->Pt() > 10){
258     // MET is needed from here
259     LoadBranch(fMetName);
260     Met *caloMet = fMet->At(0);
261     // Delta phi between the 2 leptons in degrees
262     double deltaPhiLeptons = MathUtils::DeltaPhi(leptons[0]->Phi(),
263     leptons[1]->Phi())* 180. / TMath::Pi();
264     hDFwdJetSel[11+VBYtype*100]->Fill(deltaPhiLeptons);
265     hDFwdJetSel[12+VBYtype*100]->Fill(caloMet->Pt());
266     hDFwdJetSel[13+VBYtype*100]->Fill(dilepton->Mass());
267     hDFwdJetSel[14+VBYtype*100]->Fill(nnOutputPair);
268     hDFwdJetSel[15+VBYtype*100]->Fill((double)CleanNoFwdJets->GetEntries());
269    
270     DiTauSystem ditau(leptons[0], leptons[1], caloMet);
271     hDFwdJetSel[16+VBYtype*100]->Fill(TMath::Max(TMath::Min(ditau.XTau1(),1.99),-0.99));
272     hDFwdJetSel[16+VBYtype*100]->Fill(TMath::Max(TMath::Min(ditau.XTau2(),1.99),-0.99));
273     hDFwdJetSel[17+VBYtype*100]->Fill(TMath::Min(ditau.VisMass(),399.99));
274     if(ditau.XTau1() > 0 && ditau.XTau2() > 0){
275     hDFwdJetSel[18+VBYtype*100]->Fill(TMath::Min(ditau.VisMass(),399.99));
276     hDFwdJetSel[19+VBYtype*100]->Fill(TMath::Min(ditau.RecoMass(),399.99));
277     }
278     if(dilepton->Mass() < 80 && caloMet->Pt() > 30){
279     hDFwdJetSel[20+VBYtype*100]->Fill(deltaPhiLeptons);
280     hDFwdJetSel[21+VBYtype*100]->Fill((double)CleanNoFwdJets->GetEntries());
281     }
282     if(fabs(dilepton->Mass()-91.18) < 20){
283     hDFwdJetSel[22+VBYtype*100]->Fill(caloMet->Pt());
284     hDFwdJetSel[23+VBYtype*100]->Fill((double)CleanNoFwdJets->GetEntries());
285     }
286     }
287 ceballos 1.2 delete dilepton;
288 ceballos 1.1 } // At least 2 identifed leptons
289     } // jmax >= 0
290     } // FillHist
291     } // Minimum preselection
292     else {
293     if(fFillHist == true){
294     if(fUseHighestPtJets == true){
295     if(GenqqHs->GetEntries() == 2){
296     hDFwdJetSel[6+1*100]->Fill(0.0);
297     }
298     else {
299     hDFwdJetSel[6+2*100]->Fill(0.0);
300     }
301     }
302     }
303     } // Include also this to know how many events are rejected at preselection level
304    
305     //Save Objects for Other Modules to use
306     AddObjThisEvt(CleanFwdJets, fCleanFwdJetsName.Data());
307     AddObjThisEvt(CleanNoFwdJets,fCleanNoFwdJetsName.Data());
308    
309     }
310    
311     //--------------------------------------------------------------------------------------------------
312     void FwdJetEvtSelMod::SlaveBegin()
313     {
314     // Run startup code on the computer (slave) doing the actual analysis. Here,
315     // we typically initialize histograms and other analysis objects and request
316     // branches. For this module, we request a branch of the MitTree.
317    
318     ReqBranch(fMetName, fMet);
319    
320     if(fFillHist == true){
321     char sb[200];
322     for(int j=0; j<3; j++){
323     int ind = 100 * j;
324     sprintf(sb,"hDFwdJetSel_%d",ind+0); hDFwdJetSel[ind+0] = new TH1D(sb,sb,90,0.0,180.);
325     sprintf(sb,"hDFwdJetSel_%d",ind+1); hDFwdJetSel[ind+1] = new TH1D(sb,sb,100,0.0,10.);
326     sprintf(sb,"hDFwdJetSel_%d",ind+2); hDFwdJetSel[ind+2] = new TH1D(sb,sb,200,0.0,400.);
327     sprintf(sb,"hDFwdJetSel_%d",ind+3); hDFwdJetSel[ind+3] = new TH1D(sb,sb,200,0.0,400.);
328     sprintf(sb,"hDFwdJetSel_%d",ind+4); hDFwdJetSel[ind+4] = new TH1D(sb,sb,200,0.0,4000.);
329     sprintf(sb,"hDFwdJetSel_%d",ind+5); hDFwdJetSel[ind+5] = new TH1D(sb,sb,2,-1.5,1.5);
330     sprintf(sb,"hDFwdJetSel_%d",ind+6); hDFwdJetSel[ind+6] = new TH1D(sb,sb,100,0,1);
331     sprintf(sb,"hDFwdJetSel_%d",ind+7); hDFwdJetSel[ind+7] = new TH1D(sb,sb,2,-0.5,1.5);
332     sprintf(sb,"hDFwdJetSel_%d",ind+8); hDFwdJetSel[ind+8] = new TH1D(sb,sb,2,-0.5,1.5);
333     sprintf(sb,"hDFwdJetSel_%d",ind+9); hDFwdJetSel[ind+9] = new TH1D(sb,sb,10,-0.5,9.5);
334     sprintf(sb,"hDFwdJetSel_%d",ind+10); hDFwdJetSel[ind+10] = new TH1D(sb,sb,10,-0.5,9.5);
335     sprintf(sb,"hDFwdJetSel_%d",ind+11); hDFwdJetSel[ind+11] = new TH1D(sb,sb,90,0.0,180.);
336     sprintf(sb,"hDFwdJetSel_%d",ind+12); hDFwdJetSel[ind+12] = new TH1D(sb,sb,100,0.0,200.);
337     sprintf(sb,"hDFwdJetSel_%d",ind+13); hDFwdJetSel[ind+13] = new TH1D(sb,sb,200,0.0,400.);
338     sprintf(sb,"hDFwdJetSel_%d",ind+14); hDFwdJetSel[ind+14] = new TH1D(sb,sb,100,0,1);
339     sprintf(sb,"hDFwdJetSel_%d",ind+15); hDFwdJetSel[ind+15] = new TH1D(sb,sb,10,-0.5,9.5);
340     sprintf(sb,"hDFwdJetSel_%d",ind+16); hDFwdJetSel[ind+16] = new TH1D(sb,sb,150,-1.0,2.0);
341     sprintf(sb,"hDFwdJetSel_%d",ind+17); hDFwdJetSel[ind+17] = new TH1D(sb,sb,200,0.0,400.0);
342     sprintf(sb,"hDFwdJetSel_%d",ind+18); hDFwdJetSel[ind+18] = new TH1D(sb,sb,200,0.0,400.0);
343     sprintf(sb,"hDFwdJetSel_%d",ind+19); hDFwdJetSel[ind+19] = new TH1D(sb,sb,200,0.0,400.0);
344     sprintf(sb,"hDFwdJetSel_%d",ind+20); hDFwdJetSel[ind+20] = new TH1D(sb,sb,90,0.0,180.);
345     sprintf(sb,"hDFwdJetSel_%d",ind+21); hDFwdJetSel[ind+21] = new TH1D(sb,sb,10,-0.5,9.5);
346     sprintf(sb,"hDFwdJetSel_%d",ind+22); hDFwdJetSel[ind+22] = new TH1D(sb,sb,100,0.0,200.);
347     sprintf(sb,"hDFwdJetSel_%d",ind+23); hDFwdJetSel[ind+23] = new TH1D(sb,sb,10,-0.5,9.5);
348     }
349    
350     for(int i=0; i<24; i++){
351     for(int j=0; j<3; j++){
352     AddOutput(hDFwdJetSel[i+j*100]);
353     }
354     }
355     }
356     }
357    
358     //--------------------------------------------------------------------------------------------------
359     void FwdJetEvtSelMod::SlaveTerminate()
360     {
361     // Run finishing code on the computer (slave) that did the analysis
362     }
363    
364     //--------------------------------------------------------------------------------------------------
365     void FwdJetEvtSelMod::Terminate()
366     {
367     // Run finishing code on the client computer
368     }
369    
370     //------------------------------------------------------------------------------
371     double FwdJetEvtSelMod::Testmlp_qqH(double var[6]){
372     // qqH ANN
373    
374     double OUT1 ,OUT2 , OUT3 , OUT4 , OUT5;
375     double OUT6 ,OUT7 , OUT8 , OUT9 , OUT10;
376     double OUT11 ,OUT12, OUT13, OUT14, OUT15;
377     double RIN1 ,RIN2 , RIN3 , RIN4 , RIN5;
378     double RIN6 ,RIN7 , RIN8 , RIN9 , RIN10;
379     double RIN11 ,RIN12, RIN13, RIN14, RIN15;
380    
381     OUT1 = var[0]; // PTJ1
382     OUT2 = var[1]; // PTJ2
383     OUT3 = var[2]; // DELTAETA
384     OUT4 = var[3]; // DELTAPHI
385     OUT5 = var[4]; // MJJ
386     OUT6 = var[5]; // ETAMAX*ETAMIN
387    
388    
389     // layer 2
390     RIN1 = 3.401877e-01
391     +(-1.056171e-01) * OUT1
392     +(2.830992e-01) * OUT2
393     +(2.984400e-01) * OUT3
394     +(4.116473e-01) * OUT4
395     +(-3.024486e-01) * OUT5
396     +(-1.647772e-01) * OUT6;
397     RIN2 = 2.686068e-01
398     +(-2.222253e-01) * OUT1
399     +(9.109838e-02) * OUT2
400     +(4.124089e-04) * OUT3
401     +(1.299482e-01) * OUT4
402     +(-1.341908e-01) * OUT5
403     +(-6.297962e-02) * OUT6;
404     RIN3 = 4.522297e-01
405     +(4.161951e-01) * OUT1
406     +(1.357117e-01) * OUT2
407     +(2.172969e-01) * OUT3
408     +(-3.583975e-01) * OUT4
409     +(1.069689e-01) * OUT5
410     +(-4.836994e-01) * OUT6;
411     RIN4 = -2.571132e-01
412     +(-3.627684e-01) * OUT1
413     +(3.041767e-01) * OUT2
414     +(-3.433209e-01) * OUT3
415     +(-9.905560e-02) * OUT4
416     +(-3.702095e-01) * OUT5
417     +(-3.911912e-01) * OUT6;
418     RIN5 = 4.989245e-01
419     +(-2.817431e-01) * OUT1
420     +(1.292617e-02) * OUT2
421     +(3.391084e-01) * OUT3
422     +(1.126397e-01) * OUT4
423     +(-2.039686e-01) * OUT5
424     +(1.375673e-01) * OUT6;
425     RIN6 = 2.428719e-02
426     +(-6.415707e-03) * OUT1
427     +(4.727750e-01) * OUT2
428     +(-2.074832e-01) * OUT3
429     +(2.713577e-01) * OUT4
430     +(2.674546e-02) * OUT5
431     +(2.699139e-01) * OUT6;
432     RIN7 = -9.977138e-02
433     +(3.915294e-01) * OUT1
434     +(-2.166853e-01) * OUT2
435     +(-1.475417e-01) * OUT3
436     +(3.077245e-01) * OUT4
437     +(4.190265e-01) * OUT5
438     +(-4.302447e-01) * OUT6;
439     RIN8 = 4.493271e-01
440     +(2.599535e-02) * OUT1
441     +(-4.139442e-01) * OUT2
442     +(-3.077862e-01) * OUT3
443     +(1.632269e-01) * OUT4
444     +(3.902326e-01) * OUT5
445     +(-1.511071e-01) * OUT6;
446     RIN9 = -4.361163e-01
447     +(-4.419279e-01) * OUT1
448     +(-4.473301e-02) * OUT2
449     +(-4.376293e-01) * OUT3
450     +(-2.604143e-01) * OUT4
451     +(5.240741e-01) * OUT5
452     +(4.011090e-01) * OUT6;
453     RIN10 = 3.173292e-01
454     +(-2.241760e-01) * OUT1
455     +(-6.317291e-01) * OUT2
456     +(-2.830620e-01) * OUT3
457     +(1.750255e-01) * OUT4
458     +(-4.006414e-02) * OUT5
459     +(2.725541e-01) * OUT6;
460     RIN11 = 3.396690e-02
461     +(-4.485282e-01) * OUT1
462     +(-1.537250e-02) * OUT2
463     +(4.377254e-01) * OUT3
464     +(4.355614e-01) * OUT4
465     +(3.298911e-01) * OUT5
466     +(-1.948187e-01) * OUT6;
467     RIN12 = 8.802536e-02
468     +(-1.645690e-02) * OUT1
469     +(-3.593754e-02) * OUT2
470     +(-3.300631e-01) * OUT3
471     +(1.853030e-01) * OUT4
472     +(1.364193e-02) * OUT5
473     +(6.227304e-01) * OUT6;
474     RIN13 = 3.295576e-01
475     +(-4.563790e-02) * OUT1
476     +(-2.571594e-01) * OUT2
477     +(3.939400e-01) * OUT3
478     +(-1.492818e-01) * OUT4
479     +(3.167820e-01) * OUT5
480     +(4.563062e-01) * OUT6;
481     RIN14 = 1.007358e-01
482     +(-4.626892e-03) * OUT1
483     +(9.315815e-03) * OUT2
484     +(-1.653549e-01) * OUT3
485     +(6.544170e-01) * OUT4
486     +(-1.062288e-03) * OUT5
487     +(1.809237e-01) * OUT6;
488     RIN15 = 1.842185e-01
489     +(4.109720e-01) * OUT1
490     +(-1.750934e-02) * OUT2
491     +(-2.841750e-01) * OUT3
492     +(4.502524e-01) * OUT4
493     +(4.201283e-01) * OUT5
494     +(-3.523400e-01) * OUT6;
495    
496     OUT1 = Sigmoid(RIN1);
497     OUT2 = Sigmoid(RIN2);
498     OUT3 = Sigmoid(RIN3);
499     OUT4 = Sigmoid(RIN4);
500     OUT5 = Sigmoid(RIN5);
501     OUT6 = Sigmoid(RIN6);
502     OUT7 = Sigmoid(RIN7);
503     OUT8 = Sigmoid(RIN8);
504     OUT9 = Sigmoid(RIN9);
505     OUT10 = Sigmoid(RIN10);
506     OUT11 = Sigmoid(RIN11);
507     OUT12 = Sigmoid(RIN12);
508     OUT13 = Sigmoid(RIN13);
509     OUT14 = Sigmoid(RIN14);
510     OUT15 = Sigmoid(RIN15);
511    
512     // layer 3
513     RIN1 = 4.840965e-01
514     +(1.042806e-01) * OUT1
515     +(-1.559401e-01) * OUT2
516     +(2.621663e-01) * OUT3
517     +(-2.558055e-01) * OUT4
518     +(2.465964e-01) * OUT5
519     +(-5.264297e-02) * OUT6
520     +(8.963995e-02) * OUT7
521     +(-9.445861e-02) * OUT8
522     +(-2.836257e-01) * OUT9
523     +(-2.444416e-02) * OUT10
524     +(4.221516e-01) * OUT11
525     +(-1.707438e-01) * OUT12
526     +(-1.831930e-01) * OUT13
527     +(6.496654e-01) * OUT14
528     +(-2.541290e-01) * OUT15;
529     RIN2 = -3.759727e-01
530     +(3.201696e-03) * OUT1
531     +(2.814706e-01) * OUT2
532     +(4.590240e-01) * OUT3
533     +(4.426849e-01) * OUT4
534     +(2.079947e-01) * OUT5
535     +(-1.268005e-01) * OUT6
536     +(2.240432e-01) * OUT7
537     +(-1.648751e-01) * OUT8
538     +(-2.282312e-01) * OUT9
539     +(-3.116120e-01) * OUT10
540     +(1.741257e-02) * OUT11
541     +(-1.827420e-01) * OUT12
542     +(-3.576361e-01) * OUT13
543     +(1.537631e-01) * OUT14
544     +(-4.002528e-01) * OUT15;
545     RIN3 = 3.328196e-01
546     +(-3.241346e-01) * OUT1
547     +(2.649351e-01) * OUT2
548     +(-3.956927e-01) * OUT3
549     +(4.618675e-01) * OUT4
550     +(-4.378990e-01) * OUT5
551     +(4.913286e-02) * OUT6
552     +(-2.940118e-01) * OUT7
553     +(-2.419236e-01) * OUT8
554     +(3.560494e-01) * OUT9
555     +(8.887588e-02) * OUT10
556     +(1.189877e-01) * OUT11
557     +(5.520993e-01) * OUT12
558     +(1.662827e-01) * OUT13
559     +(2.030284e-01) * OUT14
560     +(-3.767420e-01) * OUT15;
561     RIN4 = -3.021794e-01
562     +(2.426857e-02) * OUT1
563     +(-4.206145e-01) * OUT2
564     +(-3.557886e-01) * OUT3
565     +(-2.914001e-01) * OUT4
566     +(-4.993384e-02) * OUT5
567     +(3.783907e-01) * OUT6
568     +(1.474936e-01) * OUT7
569     +(3.256811e-01) * OUT8
570     +(-3.321865e-01) * OUT9
571     +(-3.979840e-01) * OUT10
572     +(5.148389e-01) * OUT11
573     +(-4.002522e-01) * OUT12
574     +(4.484992e-01) * OUT13
575     +(-3.581619e-01) * OUT14
576     +(5.721040e-01) * OUT15;
577     RIN5 = -4.385875e-01
578     +(3.705718e-01) * OUT1
579     +(-4.430307e-01) * OUT2
580     +(-4.508402e-01) * OUT3
581     +(4.231011e-01) * OUT4
582     +(5.637889e-02) * OUT5
583     +(-3.123048e-01) * OUT6
584     +(-2.918703e-01) * OUT7
585     +(-6.334352e-02) * OUT8
586     +(4.725171e-01) * OUT9
587     +(3.391282e-01) * OUT10
588     +(-1.324763e-01) * OUT11
589     +(-7.793628e-02) * OUT12
590     +(8.607008e-02) * OUT13
591     +(7.129920e-02) * OUT14
592     +(2.323856e-01) * OUT15;
593    
594     OUT1 = Sigmoid(RIN1);
595     OUT2 = Sigmoid(RIN2);
596     OUT3 = Sigmoid(RIN3);
597     OUT4 = Sigmoid(RIN4);
598     OUT5 = Sigmoid(RIN5);
599    
600     // layer 4
601     double testmlp = -6.035074e+00
602     +(-6.215130e+00) * OUT1
603     +(4.345867e+00) * OUT2
604     +(5.610218e+00) * OUT3
605     +(6.488374e+00) * OUT4
606     +(5.429897e+00) * OUT5;
607    
608     testmlp = testmlp/0.98;
609     if (testmlp<=0.0) testmlp = 0.00000;
610     else if(testmlp>=1.0) testmlp = 0.99999;
611    
612     return testmlp;
613     }
614    
615     //------------------------------------------------------------------------------
616     double FwdJetEvtSelMod::Sigmoid(double x){
617     double value = 0.;
618     if(x > 37.){
619     value = 1.;
620     }
621     else if(x < -37.){
622     value = 0.;
623     }
624     else{
625     value = 1./(1.+exp(-x));
626     }
627     return value;
628     }