ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitAna/PhysicsMod/src/FwdJetEvtSelMod.cc
Revision: 1.1
Committed: Mon Oct 6 16:59:48 2008 UTC (16 years, 7 months ago) by ceballos
Content type: text/plain
Branch: MAIN
Log Message: 
3 new analysis modules

File Contents

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