ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitAna/PhysicsMod/src/FwdJetEvtSelMod.cc
Revision: 1.4
Committed: Tue Oct 14 05:15:31 2008 UTC (16 years, 7 months ago) by loizides
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +1 -1 lines
State: FILE REMOVED
Log Message: 
Moved to MitPhysics/SelMods. Macro not yet moved. Should go into MitPhysics/HWW or something like that.

File Contents

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