Report problems to ATLAS LXR Team (with time and IP address indicated)

The LXR Cross Referencer
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ general search ]
 
 
Architecture: [ linux ]
Version: [ head ] [ nightly ] [ GaudiDev ]
  Links to LXR source navigation pages for stable releases [ 12.*.* ]   [ 13.*.* ]   [ 14.*.* ]  
001 /////////////////////////////////////////////////////////////////// 
002 // WenuNp4_006954_AlpgenJimmy_Validation.cxx 
003 // Author: Osamu Jinnouchi
004 /////////////////////////////////////////////////////////////////// 
005 
006 // STL includes
007 #include <algorithm>
008 #include <math.h>
009 #include <functional>
010 
011 // Framework includes
012 #include "GaudiKernel/MsgStream.h"
013 #include "GaudiKernel/AlgFactory.h"
014 #include "GaudiKernel/IToolSvc.h"
015 
016 #include "StoreGate/StoreGateSvc.h"
017 #include "StoreGate/DataHandle.h"
018 
019 #include "HepMC/GenEvent.h"
020 #include "HepMC/GenVertex.h"
021 #include "HepMC/GenParticle.h"
022 #include "HepMC/SimpleVector.h"
023 
024 #include "HepPDT/ParticleDataTable.hh"
025 #include "HepPDT/ParticleData.hh"
026 #include "HepPDT/TableBuilder.hh"
027 
028 #include "GaudiKernel/IIncidentSvc.h"
029 #include "GaudiKernel/Incident.h"
030 #include "GaudiKernel/IPartPropSvc.h"
031 
032 #include "McParticleKernel/ITruthParticleCnvTool.h"
033 #include "McParticleKernel/ITruthParticleFilterTool.h"
034 
035 #include "McParticleEvent/TruthParticle.h"
036 #include "McParticleEvent/TruthParticleContainer.h"
037 
038 // need for McEventCollections (OJ)
039 #include "GeneratorObjects/McEventCollection.h"
040 #include "TruthHelper/IsGenStable.h"
041 #include "TruthHelper/IsGenNonInteracting.h"
042 
043 #include "MissingETEvent/MissingEtTruth.h"
044 
045 #include "GeneratorsRTT/WenuNp4_006954_AlpgenJimmy_Validation.h"
046 
047 #include "EventInfo/EventInfo.h"
048 #include "EventInfo/EventID.h"
049 #include "EventInfo/EventType.h"
050 
051 #include "KtJet/KtUtil.h"
052 #include "KtJet/KtEvent.h"
053 #include "KtJet/KtLorentzVector.h"
054 #include "KtJet/KtJetTable.h"
055 #include "KtJet/KtDistance.h"
056 #include "KtJet/KtDistanceInterface.h"
057 #include "KtJet/KtRecom.h"
058 #include "KtJet/KtRecomInterface.h"
059 
060 //////////////////////////////////////////////////////////////////////////////////////
061 /// Constructor
062 
063 WenuNp4_006954_AlpgenJimmy_Validation::WenuNp4_006954_AlpgenJimmy_Validation(const std::string& name,
064                                        ISvcLocator* pSvcLocator) : 
065   Algorithm(name, pSvcLocator),
066   m_storeGate( "StoreGateSvc", name ),
067   m_cnvTool( ),
068   m_filterTool(0)
069 {
070 
071   const HepPDT::ParticleDataTable m_particleTable("PDG Table");
072 
073   /// switches to control the analysis through job options
074 
075   declareProperty("TruthParticles", m_truthParticlesName = "GEN_EVENT");
076   declareProperty("ConvertFromESD", m_fromESD = false );
077   declareProperty("DoAODLikeSlimming", m_AODslim = false );
078 
079   declareProperty( "CnvTool",
080                    m_cnvTool = CnvTool_t( "TruthParticleCnvTool/CnvTool", 
081                                           this ),
082                    "Handle to the tool converting a McEventCollection into a "
083                    "TruthParticleContainer" );
084   declareProperty( "FilterTool",
085                    m_filterTool = FilterTool_t( "EtaPtFilterTool", 
086                                                 this ),
087                    "AOD like filter");
088 
089 }
090 
091 
092 
093 /////////////////////////////////////////////////////////////////////////////////////
094 /// Destructor - check up memory allocation
095 /// delete any memory allocation on the heap
096 
097 WenuNp4_006954_AlpgenJimmy_Validation::~WenuNp4_006954_AlpgenJimmy_Validation() {}
098 
099 ////////////////////////////////////////////////////////////////////////////////////
100 /// Initialize
101 /// initialize StoreGate
102 
103 StatusCode WenuNp4_006954_AlpgenJimmy_Validation::initialize() 
104 {
105 
106   MsgStream msg( messageService(), name() );
107 
108   msg << MSG::INFO
109       << "Initializing WenuNp4_006954_AlpgenJimmy_Validation"
110       << endreq;
111 
112   // retrieve the storegate service
113   if ( !m_storeGate.retrieve().isSuccess() ) {
114      msg << MSG::ERROR
115          << "Unable to retrieve pointer to StoreGateSvc"
116          << endreq;
117      return StatusCode::FAILURE;
118   }
119 
120   // retrieve the truthparticle converter tool
121   if ( !m_cnvTool.retrieve().isSuccess() ) {
122     msg << MSG::ERROR
123         << "Could not retrieve the truth particle converter tool !!"
124         << endreq;
125     return StatusCode::FAILURE;
126   }
127 
128   /// Return a pointer to THistSvc
129   StatusCode sc = service("THistSvc", m_thistSvc);
130   if(sc.isFailure() ){
131     msg << MSG::ERROR
132         << "Unable to retrieve pointer to THistSvc"
133         << endreq;
134     return sc;
135   }
136 
137 
138   // Particle properties stuff
139   IPartPropSvc* p_PartPropSvc = 0;
140   static const bool CREATEIFNOTTHERE(true);
141   StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
142   if(!PartPropStatus.isSuccess() ||0==p_PartPropSvc){
143     msg<< MSG::ERROR << "Could not initialize Particle Properties service" << endreq;
144     return PartPropStatus;
145   }
146 
147   // Obtain PDT 
148   m_particleTable=p_PartPropSvc->PDT();
149 
150   /// the histograms
151   h_McElectronsN   = new TH1D("h_McElectronsN"," MC N electrons (pT>10GeV) ",5,-0.5, 4.5);
152   h_McElectronsPt  = new TH1D("h_McElectronsPt"," MC electrons pT (GeV)",50, 0.0, 200.);
153   h_McElectronsEta = new TH1D("h_McElectronsEta"," MC electrons Eta (pT>10GeV)",40, -5.0, 5.0);
154   h_McElectronsPhi = new TH1D("h_McElectronsPhi"," MC electrons Phi (pT>10GeV)",40, -3.14159, 3.14159);
155 
156   h_McMuonsN   = new TH1D("h_McMuonsN"," MC N muons (pT>10GeV) ",5,-0.5, 4.5);
157   h_McMuonsPt  = new TH1D("h_McMuonsPt"," MC muons pT (GeV)",50, 0.0, 200.);
158   h_McMuonsEta = new TH1D("h_McMuonsEta"," MC muons Eta (pT>10GeV)",40, -5.0, 5.0);
159   h_McMuonsPhi = new TH1D("h_McMuonsPhi"," MC muons Phi (pT>10GeV)",40, -3.14159, 3.14159);
160 
161   h_McTausN   = new TH1D("h_McTausN"," MC N taus (pT>10GeV) ",5,-0.5, 4.5);
162   h_McTausPt  = new TH1D("h_McTausPt"," MC taus pT (GeV)",50, 0.0, 200.);
163   h_McTausEta = new TH1D("h_McTausEta"," MC taus Eta (pT>10GeV)",40, -5.0, 5.0);
164   h_McTausPhi = new TH1D("h_McTausPhi"," MC taus Phi (pT>10GeV)",40, -3.14159, 3.14159);
165 
166   h_McPhotonsN   = new TH1D("h_McPhotonsN"," MC N photons (pT>10GeV) ",5,-0.5, 4.5);
167   h_McPhotonsPt  = new TH1D("h_McPhotonsPt"," MC photons PT (GeV)",50, 0.0, 200.);
168   h_McPhotonsEta = new TH1D("h_McPhotonsEta"," MC photons Eta (pT>10GeV)",40, -5.0, 5.0);
169   h_McPhotonsPhi = new TH1D("h_McPhotonsPhi"," MC photons Phi (pT>10GeV)",40, -3.14159, 3.14159);
170 
171   h_McBottomsN   = new TH1D("h_McBottomsN"," MC N bottoms (pT>10GeV) ",5,-0.5, 4.5);
172   h_McBottomsPt  = new TH1D("h_McBottomsPt"," MC bottoms pT (GeV)",50, 0.0, 200.);
173   h_McBottomsEta = new TH1D("h_McBottomsEta"," MC bottoms Eta (pT>10GeV)",40, -5.0, 5.0);
174   h_McBottomsPhi = new TH1D("h_McBottomsPhi"," MC bottoms Phi (pT>10GeV)",40, -3.14159, 3.14159);
175 
176   // Vector/Scalar Bosons
177 
178   //h_McZPt  = new TH1D("h_McZPt" ," MC Z Pt ",800, 0.0, 800.);
179   //h_McZy   = new TH1D("h_McZy"  ," MC Z rapidity (10GeV)",40, -5.0, 5.0);
180   //h_McZEta = new TH1D("h_McZEta"," MC Z eta (10GeV)",40, -5.0, 5.0);
181   //h_McZPhi = new TH1D("h_McZPhi"," MC Z Phi (10GeV)",40, -3.14159, 3.14159);
182 
183   //h_McRecZPt  = new TH1D("h_McRecZPt" ," MC Rec Z Pt ",800, 0.0, 800.);
184   //h_McRecZy   = new TH1D("h_McRecZy"  ," MC Rec Z rapidity (10GeV)",40, -5.0, 5.0);
185   //h_McRecZEta = new TH1D("h_McRecZEta"," MC Rec Z eta (10GeV)",40, -5.0, 5.0);
186   //h_McRecZPhi = new TH1D("h_McRecZPhi"," MC Rec Z Phi (10GeV)",40, -3.14159, 3.14159);
187 
188   h_McWPt  = new TH1D("h_McWPt" ," MC W pT (GeV)",50, 0.0, 600.);
189   h_McWy   = new TH1D("h_McWy"  ," MC W rapidity (pT>10GeV)",40, -5.0, 5.0);
190   h_McWEta = new TH1D("h_McWEta"," MC W eta (pT>10GeV)",40, -5.0, 5.0);
191   h_McWPhi = new TH1D("h_McWPhi"," MC W Phi (pT>10GeV)",40, -3.14159, 3.14159);
192 
193   //h_McHPt  = new TH1D("h_McHPt" ," MC Higgs Pt ",800, 0.0, 800.);
194   //h_McHy   = new TH1D("h_McHy"  ," MC Higgs rapidity (10GeV)",40, -5.0, 5.0);
195   //h_McHEta = new TH1D("h_McHEta"," MC Higgs eta (10GeV)",40, -5.0, 5.0);
196   //h_McHPhi = new TH1D("h_McHPhi"," MC Higgs Phi (10GeV)",40, -3.14159, 3.14159);
197 
198   // General info
199 
200   h_pdgId = new TH1D("h_pdgId","PDG ID (<50)", 100, -50, 50.);
201   h_pdgIdmid = new TH1D("h_pdgIdmid","PDG ID (<1000)", 100, -1000, 1000.);
202   h_pdgIdwid = new TH1D("h_pdgIdwid","PDG ID (<5000)", 100, -5000, 5000.);
203   h_barcode = new TH1D("h_barcode","Barcode", 1000, 0., 10000.);
204 
205   h_stablep_n = new TH1D("h_stablep_n", "Number of stable particles/event", 100, 0., 2000.);
206   h_stablec_n = new TH1D("h_stablec_n", "Number of stable charged particles/event", 100, 0., 2000.);
207   h_stablep_pt = new TH1D("h_stablep_pt", "pT (GeV) of stable particles", 150, 0., 30.);
208   h_stablec_pt = new TH1D("h_stablec_pt", "pT (GeV) of stable charged particles", 150, 0., 30.);
209 
210   h_stablep_totpx = new TH1D("h_stablep_totpx", "sum pX (GeV) of stalbe particles / event", 80, -10., 10.); 
211   h_stablep_totpy = new TH1D("h_stablep_totpy", "sum pY (GeV) of stalbe particles / event", 80, -10., 10.); 
212   h_stablep_totpz = new TH1D("h_stablep_totpz", "sum pZ (GeV) of stalbe particles / event", 80, -100., 100.); 
213 
214   // matching checks
215   h_jets_deltaR12 = new TH1D("h_jets_deltaR12", " dR of jets 1-2", 25, 0.0, 5.0);
216   h_jets_deltaR23 = new TH1D("h_jets_deltaR23", " dR of jets 2-3", 25, 0.0, 5.0);
217   h_jets_deltaR34 = new TH1D("h_jets_deltaR34", " dR of jets 3-4", 25, 0.0, 5.0);
218   
219   h_jets_d1 = new TH1D("h_jets_d1", " d_1 spectrum", 25, 0.5, 2.5);
220   h_jets_d2 = new TH1D("h_jets_d2", " d_2 spectrum", 25, 0.5, 2.5);
221   h_jets_d3 = new TH1D("h_jets_d3", " d_3 spectrum", 25, 0.5, 2.5);
222   h_jets_d4 = new TH1D("h_jets_d4", " d_4 spectrum", 25, 0.5, 2.5);
223   h_jets_d5 = new TH1D("h_jets_d5", " d_5 spectrum", 25, 0.5, 2.5);
224 
225   // Missing ET
226 
227   h_MET_Int = new TH1D("h_MET_Int","Truth MET (GeV) (Int |eta|<5)", 50, 0., 1000.);
228   h_SumET_Int = new TH1D("h_SumET_Int","Truth SumET (GeV) (Int |eta|<5)", 100, 0., 3000.);
229 
230   h_MET_NonInt = new TH1D("h_MET_NonInt","Truth MET (GeV) (NonInt all eta)", 50, 0., 1000.);
231   h_SumET_NonInt = new TH1D("h_SumET_NonInt","Truth SumET (GeV) (NonInt all eta)", 100, 0., 2000.);
232 
233   h_MET_IntCentral = new TH1D("h_MET_IntCentral","Truth MET (GeV) (IntCentral |eta|<3)", 50, 0., 1000.);
234   h_SumET_IntCentral = new TH1D("h_SumET_IntCentral","Truth SumET (GeV) (IntCentral |eta|<3)", 100, 0., 3000.);
235 
236   h_MET_IntFwd = new TH1D("h_MET_IntFwd","Truth MET (GeV) (IntFwd 3<|eta|<5)", 50, 0., 500.);
237   h_SumET_IntFwd = new TH1D("h_SumET_IntFwd","Truth SumET (GeV) (IntFwd 3<|eta|<5)", 100, 0., 1000.);
238 
239   h_MET_IntOutCover = new TH1D("h_MET_IntOutCover","Truth MET (GeV) (IntOutCover |eta|>5)", 50, 0., 200.);
240   h_SumET_IntOutCover = new TH1D("h_SumET_IntOutCover","Truth SumET (GeV) (IntOutCover |eta|>5)", 50, 0., 500.);
241   
242   h_MET_Muons= new TH1D("h_MET_Muons","Truth MET (GeV) (Muons components)", 50, 0., 500.);
243   h_SumET_Muons = new TH1D("h_SumET_Muons","Truth SumET (GeV) (Muons components)", 50, 0., 500.);
244 
245   char hname[100], htitle[100];
246 
247   h_McC4JetN   = new TH1D("h_McC4JetN"," MC N C4Jet (pT>20GeV) ",20,-0.5, 19.5);
248   for (int i=0; i<5; ++i) {
249     char m_message_pt[50];
250     char m_message[50];
251     if (i==0) {
252       sprintf(m_message_pt, "inclusive all jets");
253       sprintf(m_message, "inclusive all jets (pT>20GeV)");
254     } else {
255       sprintf(m_message_pt, "%d-th jet", i);
256       sprintf(m_message, "%d-th jet", i);
257     }
258     sprintf(hname, "h_McC4JetPt_%d",i);
259     sprintf(htitle, "MC C4Jet pT (GeV) %s", m_message_pt);
260     h_McC4JetPt[i]  = new TH1D(hname, htitle, 100, 0.0, 500.);
261 
262     sprintf(hname, "h_McC4JetEta_%d",i);
263     sprintf(htitle, "MC C4Jet Eta %s", m_message);
264     h_McC4JetEta[i] = new TH1D(hname, htitle, 40, -6.0, 6.0);
265 
266     sprintf(hname, "h_McC4JetPhi_%d",i);
267     sprintf(htitle, "MC C4Jet Phi %s", m_message);
268     h_McC4JetPhi[i] = new TH1D(hname, htitle, 40, -3.14159, 3.14159);
269   }
270 
271 
272   h_McKt4JetN   = new TH1D("h_McKt4JetN"," MC N Kt4Jet (20GeV) ",20,-0.5, 19.5);
273   for (int i=0; i<5; ++i) {
274     char m_message_pt[50];
275     char m_message[50];
276     if (i==0) {
277       sprintf(m_message_pt, "inclusive all jets");
278       sprintf(m_message, "inclusive all jets (pT>20GeV)");
279     } else {
280       sprintf(m_message_pt, "%d-th jet", i);
281       sprintf(m_message, "%d-th jet", i);
282     }
283     sprintf(hname, "h_McKt4JetPt_%d",i);
284     sprintf(htitle, "MC Kt4Jet pT (GeV) %s", m_message_pt);
285     h_McKt4JetPt[i]  = new TH1D(hname, htitle, 100, 0.0, 500.);
286 
287     sprintf(hname, "h_McKt4JetEta_%d",i);
288     sprintf(htitle, "MC Kt4Jet Eta %s", m_message);
289     h_McKt4JetEta[i] = new TH1D(hname, htitle, 40, -6.0, 6.0);
290 
291     sprintf(hname, "h_McKt4JetPhi_%d",i);
292     sprintf(htitle, "MC Kt4Jet Phi %s", m_message);
293     h_McKt4JetPhi[i] = new TH1D(hname, htitle, 40, -3.15149, 3.14159);
294   }
295 
296   sc = m_thistSvc->regHist("/AANT/MC/h_McElectronsN", h_McElectronsN);
297   sc = m_thistSvc->regHist("/AANT/MC/h_McElectronsPt", h_McElectronsPt);
298   sc = m_thistSvc->regHist("/AANT/MC/h_McElectronsEta", h_McElectronsEta);
299   sc = m_thistSvc->regHist("/AANT/MC/h_McElectronsPhi", h_McElectronsPhi);
300  
301   sc = m_thistSvc->regHist("/AANT/MC/h_McMuonsN", h_McMuonsN);
302   sc = m_thistSvc->regHist("/AANT/MC/h_McMuonsPt", h_McMuonsPt);
303   sc = m_thistSvc->regHist("/AANT/MC/h_McMuonsEta", h_McMuonsEta);
304   sc = m_thistSvc->regHist("/AANT/MC/h_McMuonsPhi", h_McMuonsPhi);
305  
306   sc = m_thistSvc->regHist("/AANT/MC/h_McTausN", h_McTausN);
307   sc = m_thistSvc->regHist("/AANT/MC/h_McTausPt", h_McTausPt);
308   sc = m_thistSvc->regHist("/AANT/MC/h_McTausEta", h_McTausEta);
309   sc = m_thistSvc->regHist("/AANT/MC/h_McTausPhi", h_McTausPhi);
310  
311   sc = m_thistSvc->regHist("/AANT/MC/h_McPhotonsN", h_McPhotonsN);
312   sc = m_thistSvc->regHist("/AANT/MC/h_McPhotonsPt", h_McPhotonsPt);
313   sc = m_thistSvc->regHist("/AANT/MC/h_McPhotonsEta", h_McPhotonsEta);
314   sc = m_thistSvc->regHist("/AANT/MC/h_McPhotonsPhi", h_McPhotonsPhi);
315 
316   sc = m_thistSvc->regHist("/AANT/MC/h_McBottomsN", h_McBottomsN);
317   sc = m_thistSvc->regHist("/AANT/MC/h_McBottomsPt", h_McBottomsPt);
318   sc = m_thistSvc->regHist("/AANT/MC/h_McBottomsEta", h_McBottomsEta);
319   sc = m_thistSvc->regHist("/AANT/MC/h_McBottomsPhi", h_McBottomsPhi);
320 
321   // sc = m_thistSvc->regHist("/AANT/MC/h_McZPt", h_McZPt);
322   // sc = m_thistSvc->regHist("/AANT/MC/h_McZy", h_McZy);
323   // sc = m_thistSvc->regHist("/AANT/MC/h_McZPhi", h_McZPhi);
324   // sc = m_thistSvc->regHist("/AANT/MC/h_McZEta", h_McZEta);
325   
326   // sc = m_thistSvc->regHist("/AANT/MC/h_McRecZPt", h_McRecZPt);
327   // sc = m_thistSvc->regHist("/AANT/MC/h_McRecZy", h_McRecZy);
328   // sc = m_thistSvc->regHist("/AANT/MC/h_McRecZPhi", h_McRecZPhi);
329   // sc = m_thistSvc->regHist("/AANT/MC/h_McRecZEta", h_McRecZEta);
330 
331   sc = m_thistSvc->regHist("/AANT/MC/h_McWPt", h_McWPt);
332   sc = m_thistSvc->regHist("/AANT/MC/h_McWy", h_McWy);
333   sc = m_thistSvc->regHist("/AANT/MC/h_McWPhi", h_McWPhi);
334   sc = m_thistSvc->regHist("/AANT/MC/h_McWEta", h_McWEta);
335 
336   // sc = m_thistSvc->regHist("/AANT/MC/h_McHPt", h_McHPt);
337   // sc = m_thistSvc->regHist("/AANT/MC/h_McHy", h_McHy);
338   // sc = m_thistSvc->regHist("/AANT/MC/h_McHPhi", h_McHPhi);
339   // sc = m_thistSvc->regHist("/AANT/MC/h_McHEta", h_McHEta);
340 
341   sc = m_thistSvc->regHist("/AANT/MC/h_pdgId", h_pdgId);
342   sc = m_thistSvc->regHist("/AANT/MC/h_pdgIdmid", h_pdgIdmid);
343   sc = m_thistSvc->regHist("/AANT/MC/h_pdgIdwid", h_pdgIdwid);
344   sc = m_thistSvc->regHist("/AANT/MC/h_barcode", h_barcode);
345 
346   sc = m_thistSvc->regHist("/AANT/MC/h_stablep_n", h_stablep_n);
347   sc = m_thistSvc->regHist("/AANT/MC/h_stablec_n", h_stablec_n);
348   sc = m_thistSvc->regHist("/AANT/MC/h_stablep_pt", h_stablep_pt);
349   sc = m_thistSvc->regHist("/AANT/MC/h_stablec_pt", h_stablec_pt);
350 
351   sc = m_thistSvc->regHist("/AANT/MC/h_stablep_totpx", h_stablep_totpx);
352   sc = m_thistSvc->regHist("/AANT/MC/h_stablep_totpy", h_stablep_totpy);
353   sc = m_thistSvc->regHist("/AANT/MC/h_stablep_totpz", h_stablep_totpz);
354 
355   sc = m_thistSvc->regHist("/AANT/MC/h_jets_deltaR12", h_jets_deltaR12);
356   sc = m_thistSvc->regHist("/AANT/MC/h_jets_deltaR23", h_jets_deltaR23);
357   sc = m_thistSvc->regHist("/AANT/MC/h_jets_deltaR34", h_jets_deltaR34);
358 
359   sc = m_thistSvc->regHist("/AANT/MC/h_jets_d1", h_jets_d1);
360   sc = m_thistSvc->regHist("/AANT/MC/h_jets_d2", h_jets_d2);
361   sc = m_thistSvc->regHist("/AANT/MC/h_jets_d3", h_jets_d3);
362   sc = m_thistSvc->regHist("/AANT/MC/h_jets_d4", h_jets_d4);
363   sc = m_thistSvc->regHist("/AANT/MC/h_jets_d5", h_jets_d5);
364  
365   sc = m_thistSvc->regHist("/AANT/MC/h_MET_Int", h_MET_Int);
366   sc = m_thistSvc->regHist("/AANT/MC/h_MET_NonInt", h_MET_NonInt);
367   sc = m_thistSvc->regHist("/AANT/MC/h_MET_IntCentral", h_MET_IntCentral);
368   sc = m_thistSvc->regHist("/AANT/MC/h_MET_IntFwd", h_MET_IntFwd);
369   sc = m_thistSvc->regHist("/AANT/MC/h_MET_IntOutCover", h_MET_IntOutCover);
370   sc = m_thistSvc->regHist("/AANT/MC/h_MET_Muons", h_MET_Muons);
371  
372   sc = m_thistSvc->regHist("/AANT/MC/h_SumET_Int", h_SumET_Int);
373   sc = m_thistSvc->regHist("/AANT/MC/h_SumET_NonInt", h_SumET_NonInt);
374   sc = m_thistSvc->regHist("/AANT/MC/h_SumET_IntCentral", h_SumET_IntCentral);
375   sc = m_thistSvc->regHist("/AANT/MC/h_SumET_IntFwd", h_SumET_IntFwd);
376   sc = m_thistSvc->regHist("/AANT/MC/h_SumET_IntOutCover", h_SumET_IntOutCover);
377   sc = m_thistSvc->regHist("/AANT/MC/h_SumET_Muons", h_SumET_Muons);
378  
379   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetN", h_McC4JetN);
380   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPt_0", h_McC4JetPt[0]);
381   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetEta_0", h_McC4JetEta[0]);
382   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPhi_0", h_McC4JetPhi[0]);
383   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPt_1", h_McC4JetPt[1]);
384   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetEta_1", h_McC4JetEta[1]);
385   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPhi_1", h_McC4JetPhi[1]);
386   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPt_2", h_McC4JetPt[2]);
387   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetEta_2", h_McC4JetEta[2]);
388   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPhi_2", h_McC4JetPhi[2]);
389   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPt_3", h_McC4JetPt[3]);
390   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetEta_3", h_McC4JetEta[3]);
391   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPhi_3", h_McC4JetPhi[3]);
392   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPt_4", h_McC4JetPt[4]);
393   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetEta_4", h_McC4JetEta[4]);
394   sc = m_thistSvc->regHist("/AANT/MC/h_McC4JetPhi_4", h_McC4JetPhi[4]);
395  
396   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetN", h_McKt4JetN);
397   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPt_0", h_McKt4JetPt[0]);
398   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetEta_0", h_McKt4JetEta[0]);
399   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPhi_0", h_McKt4JetPhi[0]);
400   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPt_1", h_McKt4JetPt[1]);
401   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetEta_1", h_McKt4JetEta[1]);
402   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPhi_1", h_McKt4JetPhi[1]);
403   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPt_2", h_McKt4JetPt[2]);
404   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetEta_2", h_McKt4JetEta[2]);
405   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPhi_2", h_McKt4JetPhi[2]);
406   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPt_3", h_McKt4JetPt[3]);
407   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetEta_3", h_McKt4JetEta[3]);
408   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPhi_3", h_McKt4JetPhi[3]);
409   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPt_4", h_McKt4JetPt[4]);
410   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetEta_4", h_McKt4JetEta[4]);
411   sc = m_thistSvc->regHist("/AANT/MC/h_McKt4JetPhi_4", h_McKt4JetPhi[4]);
412 
413 
414 
415   // OJ filter
416 
417   setupFilterTool( m_AODslim );
418 
419   // Electron-Jet Matching 
420   
421   m_elecEta = 0.;
422   m_elecPhi = 0.;
423   m_elecPt  = 0.;
424   
425 
426   return StatusCode::SUCCESS;
427 }                
428 
429 ///////////////////////////////////////////////////////////////////////////////////
430 /// Finalize - delete any memory allocation from the heap
431 
432 StatusCode WenuNp4_006954_AlpgenJimmy_Validation::finalize() 
433 {
434   MsgStream msg( messageService(), name() );
435   
436   return StatusCode::SUCCESS;
437 
438 }
439 
440 //////////////////////////////////////////////////////////////////////////////////
441 /// Execute - on event by event
442 
443 StatusCode WenuNp4_006954_AlpgenJimmy_Validation::execute() 
444 {
445   MsgStream msg( messageService(), name() );
446   StatusCode sc = StatusCode::SUCCESS;
447 
448   msg << MSG::DEBUG << "execute()" << endreq;
449 
450   /** convert on the fly TruthEvent from the ESD to TruthParticleContainer - Ketevi */
451   if ( m_fromESD ) { 
452     if ( m_cnvTool->execute().isFailure() ) {
453       msg << MSG::WARNING 
454           << "Could not convert the McEventCollection into "
455           << "a TruthParticleContainer !"
456           << endreq;
457       return StatusCode::SUCCESS;
458     }
459   }
460 
461   // create a new McEventCollection (AOD) from an old one (ex:ESD)
462   // apply some filtering criterion (if any)
463   if ( m_AODslim && !m_filterTool->execute().isSuccess() ) {
464     msg << MSG::WARNING
465         << "Could not build the AOD McEventCollection ! [algtool="
466         << m_filterTool.type() << "]"
467         << endreq;
468     return StatusCode::SUCCESS;
469   }
470 
471 
472   /** Read the TruthParticleContainer. Either the one converted from ESD TruthEvent
473      as indicated above or the default one provided by the AOD where the conversion
474      is done from GEN_AOD. -Ketevi */
475   //const TruthParticleContainer * mcParts = 0;
476   //sc = m_storeGate->retrieve( mcParts, m_truthParticlesName );
477   //if ( sc.isFailure() || 0 == mcParts ) {
478   //  msg << MSG::ERROR
479   //    << "Could not retrieve TruthParticleContainer at : "
480   //    << m_truthParticlesName
481   //    << endreq;
482   // return StatusCode::SUCCESS;
483   //}
484 
485   // -------------------------
486   // Retreive Event Info
487   // -------------------------
488   const EventInfo * evInfo =0;
489   sc = m_storeGate->retrieve( evInfo, "McEventInfo");
490   if (sc.isFailure()) {
491     msg << MSG::ERROR << "Cannot get EventInfo" << endreq;
492     return StatusCode::FAILURE;
493   }
494   //m_Run=evInfo->event_ID()->run_number();
495   //m_Event=evInfo->event_ID()->event_number();
496 
497 
498   // ------------------------------
499   // Retreive McEventCollection
500   // ------------------------------
501 
502   // Event weight
503   // (A)
504   const DataHandle<EventInfo> eventInfo, eventInfoE;
505   sc = m_storeGate->retrieve(eventInfo,eventInfoE);
506   double evtWeightA1 = (eventInfo->event_type())->mc_event_weight();
507   ///double evtWeightA2 = (eventInfoE->event_type())->mc_event_weight();
508 
509   const McEventCollection * mcCol = 0;
510   sc = m_storeGate->retrieve( mcCol, m_truthParticlesName);
511   if (sc.isFailure() || mcCol==0) {
512     msg << MSG::ERROR << "Cannot get McEventCollection" << endreq;
513     return StatusCode::FAILURE;
514   }
515   
516   TruthHelper::IsGenNonInteracting isInvis;
517   TruthHelper::IsGenStable isStable;
518 
519   McEventCollection::const_iterator itEvt;
520   HepMC::GenEvent::particle_iterator itPar;
521 
522   HepMC::GenEvent *genEvt;
523   for (itEvt = mcCol->begin();  itEvt != mcCol->end(); ++itEvt) {
524 
525     genEvt = (*itEvt);
526 
527     /// Event Weight check
528     // (B)
529     HepMC::WeightContainer weights = genEvt->weights();
530     double evtWeightB1 = weights.front(); 
531     double evtWeightB2 = weights.back(); 
532     double evtWeightSize = weights.size(); 
533     msg << MSG::INFO << " EV WEIGHT " 
534         << " (A1): " << evtWeightA1
535       //<< " (A2): " << evtWeightA2
536         << " (B1): " << evtWeightB1
537         << " (B2): " << evtWeightB2 
538         << " (Bsize): " << evtWeightSize 
539         << endreq;
540 
541     unsigned pIndex=0;
542 
543     int N_electrons = 0;
544     int N_muons = 0;
545     int N_taus = 0;
546     int N_photons = 0;
547     int N_bottoms = 0;
548     int N_stablep = 0;
549     int N_stablec = 0;
550     double Total_px = 0;
551     double Total_py = 0;
552     double Total_pz = 0;
553     
554     for(itPar=genEvt->particles_begin(); itPar!=genEvt->particles_end(); ++itPar, ++pIndex ){
555 
556 
557       // TEST
558       //msg << MSG::ALWAYS << " Begining " 
559       //  << " px " << (*itPar)->momentum().px()
560       //  << " py " << (*itPar)->momentum().py()
561       //  << " pz " << (*itPar)->momentum().pz()
562       //  << " e  " << (*itPar)->momentum().e()
563       //  << endreq;
564     
565       //double temp_px = (*itPar)->momentum().px()*1000.;
566       //double temp_py = (*itPar)->momentum().py()*1000.;
567       //double temp_pz = (*itPar)->momentum().pz()*1000.;
568       //double temp_e  = (*itPar)->momentum().e()*1000.;
569 
570       //HepMC::FourVector fv(temp_px,temp_py,temp_pz,temp_e); 
571       //(*itPar)->set_momentum(fv);
572       
573       // msg << MSG::ALWAYS << " After " 
574       //          << " px " << (*itPar)->momentum().px()
575       //          << " py " << (*itPar)->momentum().py()
576       //          << " pz " << (*itPar)->momentum().pz()
577       //          << " e  " << (*itPar)->momentum().e()
578       //          << endreq;
579 
580       const HepMC::GenParticle * mcp = (*itPar);
581       int apdgid = abs(mcp->pdg_id());
582       int bcode = (int)mcp->barcode();
583       
584       //if (1) msg << MSG::DEBUG << " Index " << pIndex 
585       //           << " PDG-ID: " << mcp->pdg_id()
586       //                 << " status: " << (isStable.operator()(mcp))
587       //                 << " invisible: " << (isInvis.operator()(mcp)) 
588       //                 << endreq;
589       
590       h_barcode -> Fill((double)bcode);
591       
592       if ( apdgid< 50. && mcp->status()==1 && !mcp->end_vertex()) {
593         h_pdgId -> Fill((float)mcp->pdg_id());
594         h_pdgIdmid -> Fill((float)mcp->pdg_id());
595         h_pdgIdwid -> Fill((float)mcp->pdg_id());
596       }
597 
598       if (mcp->status()==1 && !mcp->end_vertex()) {
599 
600         const HepPDT::ParticleData* pData = m_particleTable->particle(HepPDT::ParticleID( apdgid ));
601         //const HepMC::ParticleData* pData = m_particleTable->find( apdgid );
602         double pCharge = 0;
603         if (pData) { pCharge = pData->charge(); }
604 
605         N_stablep ++;
606         if (pCharge!=0) N_stablec ++;
607         Total_px += mcp->momentum().px()/GeV;
608         Total_py += mcp->momentum().py()/GeV;
609         Total_pz += mcp->momentum().pz()/GeV;
610         h_stablep_pt -> Fill(mcp->momentum().perp()/GeV);
611         if (pCharge!=0) h_stablec_pt -> Fill(mcp->momentum().perp()/GeV);
612       }
613 
614 
615       // Electrons
616       if (apdgid == 11 && mcp->status()==1 && !mcp->end_vertex()) {
617         h_McElectronsPt -> Fill(mcp->momentum().perp()/GeV);
618         if (mcp->momentum().perp()/GeV > 10.0) {
619           N_electrons++;
620           h_McElectronsEta -> Fill(mcp->momentum().eta());
621           h_McElectronsPhi -> Fill(mcp->momentum().phi());
622         }
623       } 
624 
625       // Muons
626       if (apdgid == 13 && mcp->status()==1 && !mcp->end_vertex()) {
627         h_McMuonsPt -> Fill(mcp->momentum().perp()/GeV);
628         if (mcp->momentum().perp()/GeV > 10.0) {
629           N_muons++;
630           h_McMuonsEta -> Fill(mcp->momentum().eta());
631           h_McMuonsPhi -> Fill(mcp->momentum().phi());
632         }
633       } 
634 
635       // Taus
636       if (apdgid == 15) {
637         HepLorentzVector pvis(0,0,0,0);
638         ///int leps = 0;
639         
640         HepMC::GenVertex *tauDvertex = mcp->end_vertex();
641         if (!tauDvertex) {
642           msg << MSG::INFO << "non decayed taus barcode: " << mcp->barcode() << endreq;
643           continue;
644         } else {
645           msg << MSG::INFO << "found decayed taus barcode: " << mcp->barcode() << endreq;
646         } 
647 
648         HepMC::GenVertex::particle_iterator curD, firstD, lastD;
649         
650         firstD = mcp->end_vertex()->particles_begin(HepMC::children);
651         lastD  = mcp->end_vertex()->particles_end(HepMC::children);
652         
653         for (curD=firstD; curD!=lastD; ++curD) {
654           int jid = abs((*curD)->pdg_id() );
655           ////if( jid>=11 && jid<=14 ) ++leps;
656           if( jid!=12 && jid!=14 && jid!=16 ) {
657             HepLorentzVector tauV;
658             HepMC::FourVector tauV4 = (*curD)->momentum();
659             tauV.setPx(tauV4.px());
660             tauV.setPy(tauV4.py());
661             tauV.setPz(tauV4.pz());
662             tauV.setE(tauV4.e());
663 
664             pvis += tauV;
665             
666             h_McTausPt -> Fill(pvis.perp()/GeV);
667             if (pvis.perp()/GeV > 10.0) {
668               N_taus++;
669               h_McTausEta -> Fill(pvis.pseudoRapidity());
670               h_McTausPhi -> Fill(pvis.phi());
671             }
672           }
673         }
674 
675 //      for(unsigned int j=0; j< mcptr.nDecay(); ++j) {
676 //        int jid = abs(mcptr.pdgDecay(j));
677 //        if( jid>=11 && jid<=14 ) ++leps;
678 //        if( jid!=12 && jid!=14 && jid!=16 ) {
679 //          pvis += mcptr.pDecay(j);
680 
681 //          h_McTausPt -> Fill(pvis.perp()/GeV);
682 //          if (pvis.perp()/GeV > 10.0) {
683 //            N_taus++;
684 //            h_McTausEta -> Fill(pvis.pseudoRapidity());
685 //            h_McTausPhi -> Fill(pvis.phi());
686 //          }
687 //        }
688 //      }
689 
690       }
691 
692       // Photons
693       if (apdgid == 22 && mcp->status()==1 && !mcp->end_vertex()) {
694         h_McPhotonsPt -> Fill(mcp->momentum().perp()/GeV);
695         if (mcp->momentum().perp()/GeV > 10.0) { 
696           N_photons++;
697           h_McPhotonsEta -> Fill(mcp->momentum().eta());
698           h_McPhotonsPhi -> Fill(mcp->momentum().phi());
699         }
700       } 
701 
702 //       // Z
703 //       if (apdgid == 23 && ((mcp->status()==195 || mcp->status()==120) // Herwig 
704 //                         || mcp->status()==155 || mcp->status()==2)) {
705         
706 //         h_McZPt -> Fill(mcp->momentum().perp()/GeV);
707 //      if (mcp->momentum().perp()/GeV > 10.) {
708 //        double myZE = mcp->momentum().e();
709 //        double myZPz= mcp->momentum().pz();
710 //        h_McZy  -> Fill(0.5*log((myZE+myZPz)/(myZE-myZPz)));
711 //        h_McZPhi -> Fill(mcp->momentum().phi());
712 //        h_McZEta -> Fill(mcp->momentum().eta());
713 //      }
714 
715 //      // goes down to stable leptons
716 //      const HepMC::GenVertex * z_decayVtx = mcp->end_vertex();
717 //      const HepMC::GenParticle *lepPlus  =0;
718 //      const HepMC::GenParticle *lepMinus =0;
719 
720 //      HepMC::GenVertex::particles_out_const_iterator l_mcpartItr  = z_decayVtx->particles_out_const_begin();
721 //      HepMC::GenVertex::particles_out_const_iterator l_mcpartItrE = z_decayVtx->particles_out_const_end();
722 
723 //      for (; l_mcpartItr != l_mcpartItrE; ++l_mcpartItr) {
724 //        HepMC::GenParticle *lep_mcpart = (*l_mcpartItr);
725 //        const HepMC::GenVertex *l_decayVtx = (*l_mcpartItr)->end_vertex();
726 //        int tmp_pdg = lep_mcpart->pdg_id();
727           
728 //        msg << MSG::DEBUG << " lepton pdg_id() " << lep_mcpart->pdg_id() << " barcode "<< lep_mcpart->barcode() << endreq;
729 //        if (l_decayVtx) msg << MSG::DEBUG 
730 //                            << " vertex barcode: "<< l_decayVtx->barcode() 
731 //                            << " vertex Id: "<< l_decayVtx->id() 
732 //                            << " vertex size: " <<l_decayVtx->particles_out_size() << endreq;
733 
734 //        int lep_iterate = 0;
735 //        while (l_decayVtx) {
736 //          HepMC::GenVertex::particles_out_const_iterator l_nextItr  = l_decayVtx->particles_out_const_begin();
737 //          HepMC::GenVertex::particles_out_const_iterator l_nextItrE = l_decayVtx->particles_out_const_end();
738 //          for (; l_nextItr != l_nextItrE; ++l_nextItr) {
739 //            msg << MSG::DEBUG << " next Itr pdgID: "<< (*l_nextItr)->pdg_id() << endreq;
740 //            if ((*l_nextItr)->pdg_id() == tmp_pdg) {
741 //              lep_mcpart = (*l_nextItr);
742 //              lep_iterate++;
743 //              break;
744 //            }
745 //          }
746 //          if (lep_iterate==0) break; // break from finding next "myself"
747 //          // next vertex point
748 //          l_decayVtx = (*l_nextItr)->end_vertex();
749 //        }
750           
751 //        msg << MSG::INFO << " Lepton Iterations: " << lep_iterate << endreq;
752 
753 //        // lepton already identified
754 //        if (tmp_pdg < 0) {
755 //          lepPlus = lep_mcpart;
756 //        } else {
757 //          lepMinus = lep_mcpart;
758 //        }
759 
760 //      }
761 
762 //      msg << MSG::DEBUG << " Plus lepton : Px " << lepPlus->momentum().px() << endreq;
763 //      msg << MSG::DEBUG << " Plus lepton : Py " << lepPlus->momentum().py() << endreq;
764 //      msg << MSG::DEBUG << " Plus lepton : Pz " << lepPlus->momentum().pz() << endreq;
765         
766 //      msg << MSG::DEBUG << " Minus lepton : Px " << lepMinus->momentum().px() << endreq;
767 //      msg << MSG::DEBUG << " Minus lepton : Py " << lepMinus->momentum().py() << endreq;
768 //      msg << MSG::DEBUG << " Minus lepton : Pz " << lepMinus->momentum().pz() << endreq;
769 
770 //      HepMC::FourVector zmom(0);
771 //      zmom.setPx(lepPlus->momentum().px() + lepMinus->momentum().px());
772 //      zmom.setPy(lepPlus->momentum().py() + lepMinus->momentum().py());
773 //      zmom.setPz(lepPlus->momentum().pz() + lepMinus->momentum().pz());
774 //      zmom.setE(lepPlus->momentum().e() + lepMinus->momentum().e());
775 
776 //      double mll = zmom.perp()/GeV;
777 //      msg << MSG::DEBUG << " Z mass calcu from leptons (GeV) " << mll << endreq;
778 //      msg << MSG::DEBUG << " Z mass from record (GeV) " << mcp->momentum().perp()/GeV << endreq;
779 //      msg << MSG::DEBUG << " Difference (calc-record) (GeV) " << mll-mcp->momentum().perp()/GeV << endreq;
780 
781 //       }
782 
783       // Bottoms
784       if (apdgid == 5 && mcp->status()==155){
785 
786         h_McBottomsPt -> Fill(mcp->momentum().perp()/GeV);
787         if (mcp->momentum().perp()/GeV > 10.0) {
788           N_bottoms++;
789           h_McBottomsEta -> Fill(mcp->momentum().eta());
790           h_McBottomsPhi -> Fill(mcp->momentum().phi());
791         }
792       }
793 
794       // W
795       if (apdgid == 24 && (mcp->status()==195||mcp->status()==155||mcp->status()==2)) {
796 
797         h_McWPt -> Fill(mcp->momentum().perp()/GeV);
798         if (mcp->momentum().perp()/GeV > 10.) {
799           double myWE = mcp->momentum().e();
800           double myWPz= mcp->momentum().pz();
801           h_McWy  -> Fill(0.5*log((myWE+myWPz)/(myWE-myWPz)));
802           h_McWPhi -> Fill(mcp->momentum().phi());
803           h_McWEta -> Fill(mcp->momentum().eta());
804         }
805 
806         // find an electron (W is just before decay, e should be found immedidately)
807         const HepMC::GenVertex * w_decayVtx = mcp->end_vertex();
808 
809         HepMC::GenVertex::particles_out_const_iterator l_mcpartItr  = w_decayVtx->particles_out_const_begin();
810         HepMC::GenVertex::particles_out_const_iterator l_mcpartItrE = w_decayVtx->particles_out_const_end();
811         
812         for (; l_mcpartItr != l_mcpartItrE; ++l_mcpartItr) {
813           HepMC::GenParticle *lep_mcpart = (*l_mcpartItr);
814           if (fabs(lep_mcpart->pdg_id()) == 11) {
815             const HepMC::GenVertex *e_decayVtx = lep_mcpart->end_vertex();
816             while (e_decayVtx) {
817               HepMC::GenVertex::particles_out_const_iterator ll_mcpartItr  = e_decayVtx->particles_out_const_begin();
818               HepMC::GenVertex::particles_out_const_iterator ll_mcpartItrE = e_decayVtx->particles_out_const_end();
819               for (; ll_mcpartItr != ll_mcpartItrE; ++ll_mcpartItr) {
820                 if (fabs((*ll_mcpartItr)->pdg_id()) == 11) {
821                   e_decayVtx = (*ll_mcpartItr)->end_vertex();
822                   lep_mcpart = (*ll_mcpartItr);
823                 }
824               }
825             }
826             
827             m_elecEta = lep_mcpart->momentum().eta();
828             m_elecPhi = lep_mcpart->momentum().phi();
829             m_elecPt  = lep_mcpart->momentum().perp();
830 //          msg << MSG::INFO << " Electron "
831 //              << " pt "  << m_elecPt 
832 //              << " eta " << m_elecEta 
833 //              << " phi " << m_elecPhi << endreq;
834             
835             break;
836           } 
837         }
838 
839       }
840 //       // H
841 //       if (apdgid == 25 && (mcp->status()==195||mcp->status()==155||mcp->status()==2)) {
842 
843 //         h_McHPt -> Fill(mcp->momentum().perp()/GeV);
844 //      if (mcp->momentum().perp()/GeV > 10.) {
845 //        double myHE = mcp->momentum().e();
846 //        double myHPz= mcp->momentum().pz();
847 //        h_McHy  -> Fill(0.5*log((myHE+myHPz)/(myHE-myHPz)));
848 //        h_McHPhi -> Fill(mcp->momentum().phi());
849 //        h_McHEta -> Fill(mcp->momentum().eta());
850 //      }
851 //       }
852     }
853 
854     h_McElectronsN -> Fill(N_electrons);
855     h_McMuonsN -> Fill(N_muons);
856     h_McTausN -> Fill(N_taus);
857     h_McPhotonsN -> Fill(N_photons);
858     h_McBottomsN -> Fill(N_bottoms);
859 
860     h_stablep_n -> Fill(N_stablep);
861     h_stablec_n -> Fill(N_stablec);
862 
863     h_stablep_totpx -> Fill(Total_px);
864     h_stablep_totpy -> Fill(Total_py);
865     h_stablep_totpz -> Fill(Total_pz);
866 
867   }
868   
869   // -------------------------------
870   // TruthJets
871   // -------------------------------
872 
873   int N_C4jets = 0;
874   int N_Kt4jets = 0;
875 
876   /// get Cone4TruthJets
877   const JetCollection*  mcC4TES;
878   IParticleContainer* mcC4NEW  = new IParticleContainer(SG::VIEW_ELEMENTS);
879   mcC4NEW->clear();
880   sc=m_storeGate->retrieve( mcC4TES, "Cone4TruthJets");
881   if( sc.isFailure()  ||  !mcC4TES ) {
882     msg << MSG::WARNING
883         << "No Cone4TruthJets container found in TDS"
884         << endreq;
885     return StatusCode::SUCCESS;
886   }
887   
888   /// iterators over the container
889   JetCollection::const_iterator mcC4jetItr  = (*mcC4TES).begin();
890   JetCollection::const_iterator mcC4jetItrE = (*mcC4TES).end();
891   int mcC4index = 0;
892   for(; mcC4jetItr != mcC4jetItrE; ++mcC4jetItr) {
893     double ptj = (*mcC4jetItr)->pt();
894     double etaj = (*mcC4jetItr)->eta();
895     double phij = (*mcC4jetItr)->phi();
896     
897     //     msg << MSG::INFO << " jet " <<mcC4index 
898     //  << " pt " << ptj 
899     //  << " eta " << etaj 
900     //  << " phi " << phij << endreq;
901 
902     // remove overlap with electron only with dR
903     if (sqrt(pow(etaj-m_elecEta, 2)+pow(phij-m_elecPhi, 2))<0.2) { 
904       msg << MSG::INFO << " C4 Jet - elec matching found " << endreq;
905       continue;
906     }
907 
908     h_McC4JetPt[0] -> Fill(ptj/GeV);
909     
910     if(ptj > 20*GeV) {
911       N_C4jets++;
912       h_McC4JetEta[0] -> Fill(etaj);
913       h_McC4JetPhi[0] -> Fill(phij);
914     }
915 
916     if (mcC4index < 4) {
917       h_McC4JetPt[mcC4index+1] -> Fill(ptj/GeV);
918       h_McC4JetEta[mcC4index+1] -> Fill(etaj);
919       h_McC4JetPhi[mcC4index+1] -> Fill(phij);
920     }
921     mcC4index++;
922     mcC4NEW->push_back(*mcC4jetItr);  // C4 Jets after overlap removal
923   }
924 
925   h_McC4JetN -> Fill(N_C4jets);
926   msg << MSG::INFO
927       << " MC cone4 jet : " << mcC4index
928       << " Size : " << (*mcC4TES).size() 
929       << endreq;
930 
931   // fill the dR dists
932   //for(mcC4jetItr  = (*mcC4TES).begin(); mcC4jetItr != mcC4jetItrE; ++mcC4jetItr) {
933   for (unsigned int i=0; i<3; ++i) {
934     if ((*mcC4NEW).size() == i+1) break;
935     HepLorentzVector j_i = (*mcC4NEW).at(i)->hlv();
936     HepLorentzVector j_j = (*mcC4NEW).at(i+1)->hlv();
937     double deltaEta = j_i.eta()-j_j.eta(); 
938     double deltaPhi = j_i.vect().deltaPhi( j_j.vect());
939     double dR = sqrt(pow(deltaEta,2)+pow(deltaPhi,2)); 
940     if (i==0) h_jets_deltaR12 -> Fill(dR);
941     if (i==1) h_jets_deltaR23 -> Fill(dR);
942     if (i==2) h_jets_deltaR34 -> Fill(dR);
943   }
944 
945   // fill the di spectrum (pseudo using C4 Jets)
946   // generate container and map for Kt LorentzVectors
947   std::vector<KtJet::KtLorentzVector>   ktLVList;
948 
949   // copy input to Kt LorentzVectors
950   IParticleContainer::const_iterator firstJet = (*mcC4NEW).begin();
951   IParticleContainer::const_iterator lastJet  = (*mcC4NEW).end();
952   for ( ; firstJet != lastJet; firstJet++ ){
953     KtJet::KtLorentzVector ktLV((*firstJet)->hlv());
954     ktLVList.push_back(ktLV);
955   }
956 
957   typedef std::string                      flag_name;
958   typedef std::map<flag_name,unsigned int> flag_store;
959 
960   KtJet::KtEvent * ktBuilder;
961   flag_store m_beamTypeFlags;                    // flags for beam type options
962   flag_store m_distSchemeFlags;                  // flags for angular def
963   flag_store m_recombSchemeFlags;                // flags for recomb scheme
964   
965   // "PP", "DeltaR", "Pt" 
966   // Reconstruction/Jet/JetRec/src/JetKtFinderTool.cxx
967   ktBuilder = new KtJet::KtEvent(ktLVList, 4, 2, 2);
968 
969   //msg << MSG::INFO << " Q 1->0:(GeV) " << sqrt(ktBuilder->getDMerge(0))/GeV << endreq;
970   //msg << MSG::INFO << " Q 2->1:(GeV) " << sqrt(ktBuilder->getDMerge(1))/GeV << endreq;
971   //msg << MSG::INFO << " Q 3->2:(GeV) " << sqrt(ktBuilder->getDMerge(2))/GeV << endreq;
972   //msg << MSG::INFO << " Q 4->3:(GeV) " << sqrt(ktBuilder->getDMerge(3))/GeV << endreq;
973   //msg << MSG::INFO << " Q 5->4:(GeV) " << sqrt(ktBuilder->getDMerge(4))/GeV << endreq;
974 
975   h_jets_d1->Fill(log10(sqrt(ktBuilder->getDMerge(0))/GeV));
976   h_jets_d2->Fill(log10(sqrt(ktBuilder->getDMerge(1))/GeV));
977   h_jets_d3->Fill(log10(sqrt(ktBuilder->getDMerge(2))/GeV));
978   h_jets_d4->Fill(log10(sqrt(ktBuilder->getDMerge(3))/GeV));
979   h_jets_d5->Fill(log10(sqrt(ktBuilder->getDMerge(4))/GeV));
980 
981   /// get Kt4TruthJets
982   const JetCollection*  mcKt4TES;
983   ///IParticleContainer* mcKt4NEW  = new IParticleContainer(SG::VIEW_ELEMENTS);
984   ///mcKt4NEW->clear();
985   sc=m_storeGate->retrieve( mcKt4TES, "Kt4TruthJets");
986   if( sc.isFailure()  ||  !mcKt4TES ) {
987     msg << MSG::WARNING
988         << "No Kt4TruthJets container found in TDS"
989         << endreq;
990     return StatusCode::SUCCESS;
991   }
992   
993   JetCollection::const_iterator mcKt4jetItr  = (*mcKt4TES).begin();
994   JetCollection::const_iterator mcKt4jetItrE = (*mcKt4TES).end();
995   int mcKt4index = 0;
996   for(; mcKt4jetItr != mcKt4jetItrE; ++mcKt4jetItr) {
997     double ptj = (*mcKt4jetItr)->pt();
998     double etaj = (*mcKt4jetItr)->eta();
999     double phij = (*mcKt4jetItr)->phi();
1000 
1001     // remove overlap with electron only with dR
1002     if (sqrt(pow(etaj-m_elecEta, 2)+pow(phij-m_elecPhi, 2))<0.2) {
1003       msg << MSG::INFO << "  C4 Jet - elec matching found " << endreq;
1004       continue;
1005     }
1006     
1007 
1008     h_McKt4JetPt[0] -> Fill(ptj/GeV);
1009     if(ptj > 20*GeV) {
1010       N_Kt4jets++;
1011       h_McKt4JetEta[0] -> Fill(etaj);
1012       h_McKt4JetPhi[0] -> Fill(phij);
1013     }
1014 
1015     if (mcKt4index < 4) {
1016       h_McKt4JetPt[mcKt4index+1] -> Fill(ptj/GeV);
1017       h_McKt4JetEta[mcKt4index+1] -> Fill(etaj);
1018       h_McKt4JetPhi[mcKt4index+1] -> Fill(phij);
1019     }
1020     mcKt4index++;
1021     ///mcKt4NEW->push_back(*mcKt4jetItr);  // Kt4 Jets after overlap removal
1022   }
1023   h_McKt4JetN -> Fill(N_Kt4jets);
1024 
1025 
1026   // TruthMissing Et variables
1027   
1028   const MissingEtTruth *missEtTruth;   // truth EtMiss   
1029   if ( (m_storeGate->retrieve(missEtTruth, "MET_Truth")).isFailure()) {
1030     msg << MSG::ERROR
1031         << "Could not retrieve Truth MissingET"
1032         << endreq;
1033   }
1034   
1035   double m_mex[6],m_mey[6],m_sumet[6]; 
1036 
1037   MissingEtTruth::TruthIndex theTruth[6] = {
1038     MissingEtTruth::Int, 
1039     MissingEtTruth::NonInt,
1040     MissingEtTruth::IntCentral,
1041     MissingEtTruth::IntFwd,
1042     MissingEtTruth::IntOutCover,
1043     MissingEtTruth::Muons,
1044   };
1045 
1046   for (int i=0; i<6;++i) { 
1047     m_mex[i] = missEtTruth -> exTruth(theTruth[i])/GeV;
1048     m_mey[i] = missEtTruth -> eyTruth(theTruth[i])/GeV;
1049     m_sumet[i] = missEtTruth -> etSumTruth(theTruth[i])/GeV;
1050   }
1051 
1052   h_MET_Int        -> Fill(sqrt(pow(m_mex[0],2)+pow(m_mey[0],2)));
1053   h_MET_NonInt     -> Fill(sqrt(pow(m_mex[1],2)+pow(m_mey[1],2)));
1054   h_MET_IntCentral -> Fill(sqrt(pow(m_mex[2],2)+pow(m_mey[2],2)));
1055   h_MET_IntFwd     -> Fill(sqrt(pow(m_mex[3],2)+pow(m_mey[3],2)));
1056   h_MET_IntOutCover-> Fill(sqrt(pow(m_mex[4],2)+pow(m_mey[4],2)));
1057   h_MET_Muons      -> Fill(sqrt(pow(m_mex[5],2)+pow(m_mey[5],2)));
1058   
1059   h_SumET_Int        -> Fill(m_sumet[0]);
1060   h_SumET_NonInt     -> Fill(m_sumet[1]);
1061   h_SumET_IntCentral -> Fill(m_sumet[2]);
1062   h_SumET_IntFwd     -> Fill(m_sumet[3]);
1063   h_SumET_IntOutCover-> Fill(m_sumet[4]);
1064   h_SumET_Muons      -> Fill(m_sumet[5]);
1065   
1066   msg << MSG::DEBUG << "EtMiss() succeeded" << endreq;
1067 
1068 //   unsigned int iPart = 0;
1069 //   for ( TruthParticleContainer::const_iterator itr = mcParts->begin();
1070 //      itr != mcParts->end();
1071 //      ++itr, ++iPart ) {
1072 //     msg << MSG::INFO << "Part " << iPart
1073 //      << " PDG-ID: " << (*itr)->pdgId()
1074 //      << " nChildren: " << (*itr)->nDecay()
1075 //      << " status: " << (*itr)->genParticle()->status()
1076 //      << " bc: " << (*itr)->genParticle()->barcode()
1077 //      << endreq;
1078 //     for ( unsigned int iChild = 0; iChild != (*itr)->nDecay(); ++iChild ){
1079 //       const TruthParticle * child = (*itr)->child( iChild );
1080 //       if ( 0 != child ) {
1081 //      msg << MSG::INFO
1082 //          << "\tchild: " << iChild
1083 //          << "\tPDGID: " << child->pdgId()
1084 //          << " status: " << child->genParticle()->status()
1085 //          << " bc: "     << child->genParticle()->barcode()
1086 //          << " bc Parents: " << child->nParents() << " [ ";
1087 //      for ( unsigned int iMoth = 0; iMoth != child->nParents(); ++iMoth ) {
1088 //        msg << child->genMother(iMoth)->barcode() << " ";
1089 //      }
1090 //      msg << "]" << endreq;
1091 //       } else {
1092 //      msg << MSG::WARNING << "Wrong pointer to child !!" << endreq;
1093 //       }
1094 //     }//> loop over children
1095 //   }//> end loop over TruthParticles
1096 
1097   return StatusCode::SUCCESS;
1098 }
1099 
1100 void WenuNp4_006954_AlpgenJimmy_Validation::setupFilterTool( Property& /*doFiltering*/ )
1101 {
1102   MsgStream m_msg( messageService(), name() );
1103     
1104   if ( m_AODslim ) {
1105     
1106     if ( !m_filterTool.retrieve().isSuccess() ) {
1107       m_msg << MSG::ERROR
1108             << "Could not retrieve algTool ITruthParticleFilterTool ["
1109             << m_filterTool.type() << "/" 
1110         //<< m_filterTool.name() 
1111             << "] !!"
1112             << endreq;
1113       throw std::runtime_error("Could not setup FilterTool property !");
1114     } else {
1115       
1116       m_msg << MSG::INFO << "Retrieved and configured algTool ["
1117             << m_filterTool.type() << "/" 
1118         //<< m_filterTool.name() 
1119             << "]"
1120             << endreq;
1121     }
1122   } else {
1123     if ( !m_filterTool.release().isSuccess() ) {
1124       m_msg << MSG::WARNING
1125             << "Could not release algTool ITruthParticleFilterTool ["
1126             << m_filterTool.type() << "/" 
1127         //<< m_filterTool.name() 
1128             << "] !!"
1129             << endreq
1130             << "Memory won't be freed until ::finalize()..."
1131             << endreq;
1132     } else {
1133       m_msg << MSG::DEBUG << "Released algTool ["
1134             << m_filterTool.type() << "/" 
1135         //<< m_filterTool.name() 
1136             << "]"
1137             << endreq;
1138     }
1139   }
1140   
1141   return;
1142 }
[ source navigation ] [ diff markup ] [ identifier search ] [ general search ]
Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine. Valid HTML 4.01!