Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-02 04:16:32

0001 /**
0002    \class HydjetAnalyzer
0003    \brief HepMC events analyzer
0004    \version 2.1
0005    \authors Yetkin Yilmaz, Andrey Belyaev
0006 */
0007 
0008 // system include files
0009 #include <cmath>
0010 #include <fstream>
0011 #include <iostream>
0012 #include <memory>
0013 #include <string>
0014 
0015 // user include files
0016 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027 
0028 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0029 
0030 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0031 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0032 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0033 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
0034 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0035 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0036 
0037 #include "HepMC/GenEvent.h"
0038 #include "HepMC/HeavyIon.h"
0039 
0040 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0041 
0042 // root include file
0043 #include "TFile.h"
0044 #include "TNtuple.h"
0045 #include "TH1.h"
0046 
0047 #include "FWCore/Framework/interface/ConsumesCollector.h"
0048 #include "FWCore/Framework/interface/EDProducer.h"
0049 #include "FWCore/Utilities/interface/EDGetToken.h"
0050 
0051 using namespace edm;
0052 using namespace std;
0053 namespace {
0054   static const int MAXPARTICLES = 5000000;
0055   static const int ETABINS = 3;  // Fix also in branch string
0056 }  // namespace
0057 struct HydjetEvent {
0058   int event;
0059   float b;
0060   float npart;
0061   float ncoll;
0062   float nhard;
0063   float phi0;
0064   float scale;
0065   int n[ETABINS];
0066   float ptav[ETABINS];
0067   int mult;
0068   float px[MAXPARTICLES];
0069   float py[MAXPARTICLES];
0070   float pz[MAXPARTICLES];
0071   float e[MAXPARTICLES];
0072   float pseudoRapidity[MAXPARTICLES];
0073   float pt[MAXPARTICLES];
0074   float eta[MAXPARTICLES];
0075   float phi[MAXPARTICLES];
0076   int pdg[MAXPARTICLES];
0077   int chg[MAXPARTICLES];
0078 
0079   float vx;
0080   float vy;
0081   float vz;
0082   float vr;
0083 };
0084 class HydjetAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0085 public:
0086   explicit HydjetAnalyzer(const edm::ParameterSet &);
0087   ~HydjetAnalyzer();
0088 
0089 private:
0090   void beginJob() final;
0091   void analyze(const edm::Event &, const edm::EventSetup &) final;
0092   void endJob() final;
0093   // ----------member data ---------------------------
0094   std::ofstream out_b;
0095   std::string fBFileName;
0096   std::ofstream out_n;
0097   std::string fNFileName;
0098   std::ofstream out_m;
0099   std::string fMFileName;
0100   TTree *hydjetTree_;
0101   HydjetEvent hev_;
0102   TNtuple *nt;
0103   std::string output;  // Output filename
0104   bool doAnalysis_;
0105   bool printLists_;
0106   bool doCF_;
0107   bool doVertex_;
0108   bool useHepMCProduct_;
0109   bool doHI_;
0110   bool doParticles_;
0111   double etaMax_;
0112   double ptMin_;
0113   bool doHistos_, userHistos_;
0114 
0115   float *ptBins;
0116   float *etaBins;
0117   float *phiBins;
0118   float *v2ptBins;
0119   float *v2etaBins;
0120 
0121   vector<double> uPtBins_;
0122   vector<double> uEtaBins_;
0123   vector<double> uPhiBins_;
0124   vector<double> uV2ptBins_;
0125   vector<double> uV2etaBins_;
0126 
0127   int uPDG_1;
0128   int uPDG_2;
0129   int uPDG_3;
0130   int uStatus_;
0131   int nintPt = 0;
0132   int nintEta = 0;
0133   int nintPhi = 0;
0134   int nintV2pt = 0;
0135   int nintV2eta = 0;
0136 
0137   double minPt = 0.;
0138   double minEta = 0.;
0139   double minPhi = 0.;
0140   double minV2pt = 0.;
0141   double minV2eta = 0.;
0142 
0143   double maxPt = 0.;
0144   double maxEta = 0.;
0145   double maxPhi = 0.;
0146   double maxV2pt = 0.;
0147   double maxV2eta = 0.;
0148 
0149   double upTetaCut_ = 0.;
0150   double downTetaCut_ = -1.;
0151   const double pi = 3.14159265358979;
0152 
0153   edm::EDGetTokenT<edm::HepMCProduct> srcT_;
0154   edm::EDGetTokenT<CrossingFrame<edm::HepMCProduct>> srcTmix_;
0155 
0156   edm::InputTag genParticleSrc_;
0157   edm::InputTag genHIsrc_;
0158   edm::InputTag simVerticesTag_;
0159   edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> pdtToken_;
0160 
0161   //common
0162   TH1D *dhpt;
0163   TH1D *dhpt_ch;
0164 
0165   TH1D *dhphi;
0166   TH1D *dhphi_ch;
0167 
0168   TH1D *dhv2pt;
0169   TH1D *dhv0pt;
0170 
0171   TH1D *dhv2eta;
0172   TH1D *dhv0eta;
0173 
0174   TH1D *dheta;
0175   TH1D *dheta_ch;
0176 
0177   TH1D *hNev;
0178 
0179   //Users
0180   TH1D *uhpt;
0181   TH1D *uhpth;
0182   TH1D *uhptj;
0183 
0184   TH1D *uhpt_db;
0185   TH1D *uhpth_db;
0186   TH1D *uhptj_db;
0187 
0188   TH1D *uhNpart;
0189   TH1D *uhNparth;
0190   TH1D *uhNpartj;
0191 
0192   TH1D *uhNpart_db;
0193   TH1D *uhNparth_db;
0194   TH1D *uhNpartj_db;
0195 
0196   TH1D *uhPtNpart;
0197   TH1D *uhPtNparth;
0198   TH1D *uhPtNpartj;
0199 
0200   TH1D *uhPtNpart_db;
0201   TH1D *uhPtNparth_db;
0202   TH1D *uhPtNpartj_db;
0203 
0204   TH1D *uhv2Npart;
0205   TH1D *uhv2Nparth;
0206   TH1D *uhv2Npartj;
0207 
0208   TH1D *uhv2Npart_db;
0209   TH1D *uhv2Nparth_db;
0210   TH1D *uhv2Npartj_db;
0211 
0212   TH1D *uheta;
0213   TH1D *uhetah;
0214   TH1D *uhetaj;
0215 
0216   TH1D *uhphi;
0217   TH1D *uhphih;
0218   TH1D *uhphij;
0219 
0220   TH1D *uhv2pt;
0221   TH1D *uhv2pth;
0222   TH1D *uhv2ptj;
0223 
0224   TH1D *uhv3pt;
0225   TH1D *uhv4pt;
0226   TH1D *uhv5pt;
0227   TH1D *uhv6pt;
0228 
0229   TH1D *uhv0pt;
0230   TH1D *uhv0pth;
0231   TH1D *uhv0ptj;
0232 
0233   TH1D *uhv2pt_db;
0234   TH1D *uhv2pth_db;
0235   TH1D *uhv2ptj_db;
0236 
0237   TH1D *uhv3pt_db;
0238   TH1D *uhv4pt_db;
0239   TH1D *uhv5pt_db;
0240   TH1D *uhv6pt_db;
0241 
0242   TH1D *uhv0pt_db;
0243   TH1D *uhv0pth_db;
0244   TH1D *uhv0ptj_db;
0245 
0246   TH1D *uhv2eta;
0247   TH1D *uhv2etah;
0248   TH1D *uhv2etaj;
0249 
0250   TH1D *uhv3eta;
0251   TH1D *uhv4eta;
0252   TH1D *uhv5eta;
0253   TH1D *uhv6eta;
0254 
0255   TH1D *uhv0eta;
0256   TH1D *uhv0etah;
0257   TH1D *uhv0etaj;
0258 
0259   TH1D *uhv2eta_db;
0260   TH1D *uhv2etah_db;
0261   TH1D *uhv2etaj_db;
0262 
0263   TH1D *uhv3eta_db;
0264   TH1D *uhv4eta_db;
0265   TH1D *uhv5eta_db;
0266   TH1D *uhv6eta_db;
0267 
0268   TH1D *uhv0eta_db;
0269   TH1D *uhv0etah_db;
0270   TH1D *uhv0etaj_db;
0271 };
0272 
0273 HydjetAnalyzer::HydjetAnalyzer(const edm::ParameterSet &iConfig) {
0274   fBFileName = iConfig.getUntrackedParameter<std::string>("output_b", "b_values.txt");
0275   fNFileName = iConfig.getUntrackedParameter<std::string>("output_n", "n_values.txt");
0276   fMFileName = iConfig.getUntrackedParameter<std::string>("output_m", "m_values.txt");
0277   doAnalysis_ = iConfig.getUntrackedParameter<bool>("doAnalysis", false);
0278   useHepMCProduct_ = iConfig.getUntrackedParameter<bool>("useHepMCProduct", true);
0279   printLists_ = iConfig.getUntrackedParameter<bool>("printLists", false);
0280   doCF_ = iConfig.getUntrackedParameter<bool>("doMixed", false);
0281   doVertex_ = iConfig.getUntrackedParameter<bool>("doVertex", false);
0282   if (doVertex_) {
0283     simVerticesTag_ = iConfig.getParameter<edm::InputTag>("simVerticesTag");
0284   }
0285   etaMax_ = iConfig.getUntrackedParameter<double>("etaMax", 2.);
0286   ptMin_ = iConfig.getUntrackedParameter<double>("ptMin", 0);
0287   srcT_ = mayConsume<HepMCProduct>(
0288       iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("generator", "unsmeared")));
0289   srcTmix_ = consumes<CrossingFrame<edm::HepMCProduct>>(
0290       iConfig.getUntrackedParameter<edm::InputTag>("srcMix", edm::InputTag("mix", "generatorSmeared")));
0291 
0292   genParticleSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("hiGenParticles"));
0293   genHIsrc_ = iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("heavyIon"));
0294 
0295   if (useHepMCProduct_)
0296     pdtToken_ = esConsumes();
0297   doParticles_ = iConfig.getUntrackedParameter<bool>("doParticles", false);
0298   doHistos_ = iConfig.getUntrackedParameter<bool>("doHistos", false);
0299   if (doHistos_) {
0300     userHistos_ = iConfig.getUntrackedParameter<bool>("userHistos", false);
0301     if (userHistos_) {
0302       uStatus_ = iConfig.getUntrackedParameter<int>("uStatus");
0303       uPDG_1 = iConfig.getUntrackedParameter<int>("uPDG_1");
0304       uPDG_2 = iConfig.getUntrackedParameter<int>("uPDG_2", uPDG_1);
0305       uPDG_3 = iConfig.getUntrackedParameter<int>("uPDG_3", uPDG_1);
0306       upTetaCut_ = iConfig.getUntrackedParameter<double>("uPTetaCut", 0.8);
0307       downTetaCut_ = iConfig.getUntrackedParameter<double>("dPTetaCut", -1.);
0308       uPtBins_ = iConfig.getUntrackedParameter<vector<double>>("PtBins");
0309       uEtaBins_ = iConfig.getUntrackedParameter<vector<double>>("EtaBins");
0310       uPhiBins_ = iConfig.getUntrackedParameter<vector<double>>("PhiBins");
0311       uV2ptBins_ = iConfig.getUntrackedParameter<vector<double>>("v2PtBins");
0312       uV2etaBins_ = iConfig.getUntrackedParameter<vector<double>>("v2EtaBins");
0313 
0314       //Pt
0315       int PtSize = uPtBins_.size();
0316       if (PtSize > 1) {
0317         ptBins = new float[PtSize];
0318         nintPt = PtSize - 1;
0319         for (int k = 0; k < PtSize; k++) {
0320           ptBins[k] = uPtBins_[k];
0321         }
0322       } else {
0323         nintPt = iConfig.getUntrackedParameter<int>("nintPt");
0324         maxPt = iConfig.getUntrackedParameter<double>("maxPt");
0325         minPt = iConfig.getUntrackedParameter<double>("minPt");
0326         ptBins = new float[nintPt + 1];
0327         for (int k = 0; k < nintPt + 1; k++) {
0328           ptBins[k] = minPt + k * ((maxPt - minPt) / nintPt);
0329         }
0330       }
0331 
0332       //Eta
0333       int EtaSize = uEtaBins_.size();
0334       if (EtaSize > 1) {
0335         etaBins = new float[EtaSize];
0336         nintEta = EtaSize - 1;
0337         for (int k = 0; k < EtaSize; k++) {
0338           etaBins[k] = uEtaBins_[k];
0339         }
0340       } else {
0341         nintEta = iConfig.getUntrackedParameter<int>("nintEta");
0342         maxEta = iConfig.getUntrackedParameter<double>("maxEta");
0343         minEta = iConfig.getUntrackedParameter<double>("minEta");
0344         etaBins = new float[nintEta + 1];
0345         for (int k = 0; k < nintEta + 1; k++) {
0346           etaBins[k] = minEta + k * ((maxEta - minEta) / nintEta);
0347         }
0348       }
0349 
0350       //Phi
0351       int PhiSize = uPhiBins_.size();
0352       if (PhiSize > 1) {
0353         phiBins = new float[PhiSize];
0354         nintPhi = PhiSize - 1;
0355         for (int k = 0; k < PhiSize; k++) {
0356           phiBins[k] = uPhiBins_[k];
0357         }
0358       } else {
0359         nintPhi = iConfig.getUntrackedParameter<int>("nintPhi");
0360         maxPhi = iConfig.getUntrackedParameter<double>("maxPhi");
0361         minPhi = iConfig.getUntrackedParameter<double>("minPhi");
0362         phiBins = new float[nintPhi + 1];
0363         for (int k = 0; k < nintPhi + 1; k++) {
0364           phiBins[k] = minPhi + k * ((maxPhi - minPhi) / nintPhi);
0365         }
0366       }
0367 
0368       //v2Pt
0369       int v2PtSize = uV2ptBins_.size();
0370       if (v2PtSize > 1) {
0371         v2ptBins = new float[v2PtSize];
0372         nintV2pt = v2PtSize - 1;
0373         for (int k = 0; k < v2PtSize; k++) {
0374           v2ptBins[k] = uV2ptBins_[k];
0375         }
0376       } else {
0377         nintV2pt = iConfig.getUntrackedParameter<int>("nintV2pt");
0378         maxV2pt = iConfig.getUntrackedParameter<double>("maxV2pt");
0379         minV2pt = iConfig.getUntrackedParameter<double>("minV2pt");
0380         v2ptBins = new float[nintV2pt + 1];
0381         for (int k = 0; k < nintV2pt + 1; k++) {
0382           v2ptBins[k] = minV2pt + k * ((maxV2pt - minV2pt) / nintV2pt);
0383         }
0384       }
0385 
0386       //v2Eta
0387       int v2EtaSize = uV2etaBins_.size();
0388       if (v2EtaSize > 1) {
0389         v2etaBins = new float[v2EtaSize];
0390         nintV2eta = v2EtaSize - 1;
0391         for (int k = 0; k < v2EtaSize; k++) {
0392           v2etaBins[k] = uV2etaBins_[k];
0393         }
0394       } else {
0395         nintV2eta = iConfig.getUntrackedParameter<int>("nintV2eta");
0396         maxV2eta = iConfig.getUntrackedParameter<double>("maxV2eta");
0397         minV2eta = iConfig.getUntrackedParameter<double>("minV2eta");
0398         v2etaBins = new float[nintV2eta + 1];
0399         for (int k = 0; k < nintV2eta + 1; k++) {
0400           v2etaBins[k] = minV2eta + k * ((maxV2eta - minV2eta) / nintV2eta);
0401         }
0402       }
0403     }  //user histo
0404   }    //do histo
0405 }
0406 HydjetAnalyzer::~HydjetAnalyzer() {}
0407 
0408 // ------------ method called to for each event  ------------
0409 void HydjetAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0410   using namespace edm;
0411   using namespace HepMC;
0412   hev_.event = iEvent.id().event();
0413   for (int ieta = 0; ieta < ETABINS; ++ieta) {
0414     hev_.n[ieta] = 0;
0415     hev_.ptav[ieta] = 0;
0416   }
0417   hev_.mult = 0;
0418   double phi0 = 0.;
0419   double phi3 = 0.;
0420   double b = -1.;
0421   double v2, v3, v4, v5, v6;
0422   double scale = -1.;
0423   int npart = -1;
0424   int ncoll = -1;
0425   int nhard = -1;
0426   double vx = -99.;
0427   double vy = -99.;
0428   double vz = -99.;
0429   double vr = -99.;
0430   const GenEvent *evt;
0431   int nmix = -1;
0432   int np = 0;
0433   int sig = -1;
0434   int src = -1;
0435   if (useHepMCProduct_) {
0436     edm::ESHandle<ParticleDataTable> pdt = iSetup.getHandle(pdtToken_);
0437 
0438     if (doCF_) {
0439       Handle<CrossingFrame<HepMCProduct>> cf;
0440       iEvent.getByToken(srcTmix_, cf);
0441       MixCollection<HepMCProduct> mix(cf.product());
0442       nmix = mix.size();
0443       //cout << "Mix Collection Size: " << mix.size() <<", pileup size: " <<mix.sizePileup() << ", signal: "<<mix.sizeSignal()<< endl;
0444       MixCollection<HepMCProduct>::iterator mbegin = mix.begin();
0445       MixCollection<HepMCProduct>::iterator mend = mix.end();
0446       for (MixCollection<HepMCProduct>::iterator mixit = mbegin; mixit != mend; ++mixit) {
0447         const GenEvent *subevt = (*mixit).GetEvent();
0448         int all = subevt->particles_size();
0449         //cout << "Subevent size: " << all << " Subevent type (1-signal): "<< (mixit).getTrigger()<<" Source type (pileup=0, cosmics=1, beam halo+ =2, beam halo- =3): "<< (mixit).getSourceType()<<" Bunchcrossing number: "<< mixit.bunch()<<" Impact: " <<subevt->heavy_ion()->impact_parameter()<<endl;
0450         np += all;
0451         HepMC::GenEvent::particle_const_iterator begin = subevt->particles_begin();
0452         HepMC::GenEvent::particle_const_iterator end = subevt->particles_end();
0453         for (HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it) {
0454           if ((*it)->status() == 1) {
0455             int pdg_id = (*it)->pdg_id();
0456             float eta = (*it)->momentum().eta();
0457             float phi = (*it)->momentum().phi();
0458             float pt = (*it)->momentum().perp();
0459             const ParticleData *part = pdt->particle(pdg_id);
0460             int charge = static_cast<int>(part->charge());
0461             hev_.pt[hev_.mult] = pt;
0462             hev_.eta[hev_.mult] = eta;
0463             hev_.phi[hev_.mult] = phi;
0464             hev_.pdg[hev_.mult] = pdg_id;
0465             hev_.chg[hev_.mult] = charge;
0466 
0467             //cout << "Mix Particles: pt= " << pt<<" eta="<<eta<<" phi="<<phi<<" pdg="<< pdg_id<<" charge="<<charge << endl;
0468             eta = fabs(eta);
0469             int etabin = 0;
0470             if (eta > 0.5)
0471               etabin = 1;
0472             if (eta > 1.)
0473               etabin = 2;
0474             if (eta < 2.) {
0475               hev_.ptav[etabin] += pt;
0476               ++(hev_.n[etabin]);
0477             }
0478             ++(hev_.mult);
0479           }
0480         }
0481       }
0482     } else {  //not mixing
0483       Handle<HepMCProduct> mc;
0484       iEvent.getByToken(srcT_, mc);
0485       evt = mc->GetEvent();
0486       scale = evt->event_scale();
0487       const HeavyIon *hi = evt->heavy_ion();
0488 
0489       if (hi) {
0490         b = hi->impact_parameter();
0491         npart = hi->Npart_proj() + hi->Npart_targ();
0492         ncoll = hi->Ncoll();
0493         nhard = hi->Ncoll_hard();
0494         phi0 = hi->event_plane_angle();
0495         phi3 =
0496             hi->eccentricity();  // 0.;  // No HepMC entry for Psi3 exist, but in private code it's possible to use hi->eccentricity();
0497         if (printLists_) {
0498           out_b << b << endl;
0499           out_n << npart << endl;
0500         }
0501       }
0502 
0503       src = evt->particles_size();
0504 
0505       HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
0506       HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
0507       for (HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it) {
0508         if (((*it)->status() > 0) && ((*it)->status() < 21)) {
0509           const ParticleData *part;
0510           int pdg_id = (*it)->pdg_id();
0511           float eta = (*it)->momentum().eta();
0512           float phi = (*it)->momentum().phi();
0513           float pt = (*it)->momentum().perp();
0514 
0515           float px = (*it)->momentum().px();
0516           float py = (*it)->momentum().py();
0517           float pz = (*it)->momentum().pz();
0518           float e = (*it)->momentum().e();
0519           float pseudoRapidity = (*it)->momentum().pseudoRapidity();
0520           int charge = -1;
0521           if ((pdt->particle(pdg_id))) {
0522             part = pdt->particle(pdg_id);
0523             charge = static_cast<int>(part->charge());
0524           }
0525 
0526           /*
0527           if(pdg_id==-130){ //there are not -130 in pdt
0528         part = pdt.particle(130);
0529       }else if(!(pdt.particle(pdg_id))){ //skip if not in PDT!!!
0530         cout<<" Out of PDT: "<< pdg_id<<endl;
0531         continue;
0532           }else{
0533         part = pdt.particle(pdg_id);
0534           }
0535 */
0536           //          int charge = static_cast<int>(part->charge());
0537           hev_.px[hev_.mult] = px;
0538           hev_.py[hev_.mult] = py;
0539           hev_.pz[hev_.mult] = pz;
0540           hev_.e[hev_.mult] = e;
0541           hev_.pseudoRapidity[hev_.mult] = pseudoRapidity;
0542           hev_.pt[hev_.mult] = pt;
0543           hev_.eta[hev_.mult] = eta;
0544           hev_.phi[hev_.mult] = phi;
0545           hev_.pdg[hev_.mult] = pdg_id;
0546           hev_.chg[hev_.mult] = charge;
0547 
0548           phi = phi - phi0;
0549           ///
0550           double phiTrue;
0551           if (phi > pi) {
0552             phiTrue = phi - (2 * pi);
0553           } else if (phi < (-1 * pi)) {
0554             phiTrue = phi + (2 * pi);
0555           } else {
0556             phiTrue = phi;
0557           }
0558           ///
0559           v2 = std::cos(2 * (phiTrue));
0560           v3 = std::cos(3 * (phiTrue - phi3));
0561           v4 = std::cos(4 * (phiTrue));
0562           v5 = std::cos(5 * (phiTrue - phi3));
0563           v6 = std::cos(6 * (phiTrue));
0564 
0565           if (doHistos_) {
0566             //common histos
0567             if ((*it)->status() < 10) {  //status 1
0568               dheta->Fill(eta);
0569 
0570               if (std::abs(eta) < 0.8) {
0571                 dhpt->Fill(pt);
0572                 dhphi->Fill(phiTrue);
0573               }
0574 
0575               if (std::abs(pdg_id) == 211 || std::abs(pdg_id) == 321 || std::abs(pdg_id) == 2212) {  //ch
0576 
0577                 if (std::abs(eta) < 0.8) {
0578                   dhv0pt->Fill(pt, 1.);
0579                   dhv2pt->Fill(pt, v2);
0580                   dhpt_ch->Fill(pt);
0581                   dhphi_ch->Fill(phiTrue);
0582                 }
0583 
0584                 dhv0eta->Fill(eta, 1.);
0585                 dhv2eta->Fill(eta, v2);
0586                 dheta_ch->Fill(eta);
0587               }  //ch
0588             }    //status 1
0589 
0590             //user histos
0591             if (userHistos_ && ((uStatus_ == 3) || (((*it)->status() < 10) && (uStatus_ == 1)) ||
0592                                 (((*it)->status() > 10) && (uStatus_ == 2)))) {  //user status
0593 
0594               //set1
0595               if (std::abs(pdg_id) == uPDG_1 || std::abs(pdg_id) == uPDG_2 || std::abs(pdg_id) == uPDG_3) {  //uPDG
0596                 if ((uStatus_ == 3) && ((*it)->status() < 10))
0597                   cout << "ustatus=3, but stab. part. found!!!" << endl;
0598 
0599                 if (std::abs(eta) > downTetaCut_ && std::abs(eta) < upTetaCut_) {  //eta cut
0600 
0601                   uhv0pt->Fill(pt, 1.);
0602                   uhv2pt->Fill(pt, v2);
0603                   uhv3pt->Fill(pt, v3);
0604                   uhv4pt->Fill(pt, v4);
0605                   uhv5pt->Fill(pt, v5);
0606                   uhv6pt->Fill(pt, v6);
0607 
0608                   uhv0pt_db->Fill(pt, 1.);
0609                   uhv2pt_db->Fill(pt, v2);
0610                   uhv3pt_db->Fill(pt, v3);
0611                   uhv4pt_db->Fill(pt, v4);
0612                   uhv5pt_db->Fill(pt, v5);
0613                   uhv6pt_db->Fill(pt, v6);
0614 
0615                   if (pt >= 1.5 && pt < 10.) {
0616                     uhv2Npart->Fill(npart, v2);
0617                     uhv2Npart_db->Fill(npart, v2);
0618 
0619                     uhPtNpart->Fill(npart, pt);
0620                     uhPtNpart_db->Fill(npart, pt);
0621 
0622                     uhNpart->Fill(npart, 1.);
0623                     uhNpart_db->Fill(npart, 1.);
0624                   }
0625 
0626                   uhpt->Fill(pt);
0627                   uhpt_db->Fill(pt);
0628                   uhphi->Fill(phiTrue);
0629 
0630                   if (((*it)->status() == 16) || ((*it)->status() == 6)) {  //hydro
0631                     uhv0pth->Fill(pt, 1.);
0632                     uhv2pth->Fill(pt, v2);
0633 
0634                     uhv0pth_db->Fill(pt, 1.);
0635                     uhv2pth_db->Fill(pt, v2);
0636 
0637                     if (pt >= 1.5 && pt < 10.) {
0638                       uhv2Nparth->Fill(npart, v2);
0639                       uhv2Nparth_db->Fill(npart, v2);
0640                     }
0641 
0642                     uhPtNparth->Fill(npart, pt);
0643                     uhPtNparth_db->Fill(npart, pt);
0644 
0645                     uhpth->Fill(pt);
0646                     uhpth_db->Fill(pt);
0647                     uhphih->Fill(phiTrue);
0648                   }
0649 
0650                   if (((*it)->status() == 17) || ((*it)->status() == 7)) {  //jet
0651                     uhv0ptj->Fill(pt, 1.);
0652                     uhv2ptj->Fill(pt, v2);
0653 
0654                     uhv0ptj_db->Fill(pt, 1.);
0655                     uhv2ptj_db->Fill(pt, v2);
0656 
0657                     if (pt >= 1.5 && pt < 10.) {
0658                       uhv2Npartj->Fill(npart, v2);
0659                       uhv2Npartj_db->Fill(npart, v2);
0660                     }
0661 
0662                     uhPtNpartj->Fill(npart, pt);
0663                     uhPtNpartj_db->Fill(npart, pt);
0664 
0665                     uhptj->Fill(pt);
0666                     uhptj_db->Fill(pt);
0667                     uhphij->Fill(phiTrue);
0668                   }
0669                 }  //eta cut
0670 
0671                 uheta->Fill(eta);
0672 
0673                 uhv0eta->Fill(eta, 1.);
0674                 uhv2eta->Fill(eta, v2);
0675                 uhv3eta->Fill(eta, v3);
0676                 uhv4eta->Fill(eta, v4);
0677                 uhv5eta->Fill(eta, v5);
0678                 uhv6eta->Fill(eta, v6);
0679 
0680                 uhv0eta_db->Fill(eta, 1.);
0681                 uhv2eta_db->Fill(eta, v2);
0682                 uhv3eta_db->Fill(eta, v3);
0683                 uhv4eta_db->Fill(eta, v4);
0684                 uhv5eta_db->Fill(eta, v5);
0685                 uhv6eta_db->Fill(eta, v6);
0686 
0687                 if (((*it)->status() == 16) || ((*it)->status() == 6)) {  //hydro
0688                   uhv2etah->Fill(eta, v2);
0689                   uhv0etah->Fill(eta, 1.);
0690 
0691                   uhv2etah_db->Fill(eta, v2);
0692                   uhv0etah_db->Fill(eta, 1.);
0693 
0694                   uhetah->Fill(eta);
0695                 }
0696                 if (((*it)->status() == 17) || ((*it)->status() == 7)) {  //jet
0697                   uhv2etaj->Fill(eta, v2);
0698                   uhv0etaj->Fill(eta, 1.);
0699 
0700                   uhv2etaj_db->Fill(eta, v2);
0701                   uhv0etaj_db->Fill(eta, 1.);
0702 
0703                   uhetaj->Fill(eta);
0704                 }
0705               }  //uPDG
0706 
0707             }  //user status
0708 
0709           }  //doHistos_
0710 
0711           eta = fabs(eta);
0712           int etabin = 0;
0713           if (eta > 0.5)
0714             etabin = 1;
0715           if (eta > 1.)
0716             etabin = 2;
0717           if (eta < 2.) {
0718             hev_.ptav[etabin] += pt;
0719             ++(hev_.n[etabin]);
0720           }
0721           ++(hev_.mult);
0722         }
0723       }
0724     }
0725   } else {  // not HepMC
0726     edm::Handle<reco::GenParticleCollection> parts;
0727     iEvent.getByLabel(genParticleSrc_, parts);
0728     for (unsigned int i = 0; i < parts->size(); ++i) {
0729       const reco::GenParticle &p = (*parts)[i];
0730       hev_.pt[hev_.mult] = p.pt();
0731       hev_.eta[hev_.mult] = p.eta();
0732       hev_.phi[hev_.mult] = p.phi();
0733       hev_.pdg[hev_.mult] = p.pdgId();
0734       hev_.chg[hev_.mult] = p.charge();
0735       double eta = fabs(p.eta());
0736 
0737       int etabin = 0;
0738       if (eta > 0.5)
0739         etabin = 1;
0740       if (eta > 1.)
0741         etabin = 2;
0742       if (eta < 2.) {
0743         hev_.ptav[etabin] += p.pt();
0744         ++(hev_.n[etabin]);
0745       }
0746       ++(hev_.mult);
0747     }
0748     if (doHI_) {
0749       edm::Handle<GenHIEvent> higen;
0750       iEvent.getByLabel(genHIsrc_, higen);
0751     }
0752   }
0753 
0754   if (doVertex_) {
0755     edm::Handle<edm::SimVertexContainer> simVertices;
0756     iEvent.getByLabel<edm::SimVertexContainer>(simVerticesTag_, simVertices);
0757 
0758     if (!simVertices.isValid())
0759       throw cms::Exception("FatalError") << "No vertices found\n";
0760     int inum = 0;
0761 
0762     edm::SimVertexContainer::const_iterator it = simVertices->begin();
0763     SimVertex vertex = (*it);
0764     cout << " Vertex position " << inum << " " << vertex.position().rho() << " " << vertex.position().z() << endl;
0765     vx = vertex.position().x();
0766     vy = vertex.position().y();
0767     vz = vertex.position().z();
0768     vr = vertex.position().rho();
0769   }
0770 
0771   for (int i = 0; i < 3; ++i) {
0772     hev_.ptav[i] = hev_.ptav[i] / hev_.n[i];
0773   }
0774 
0775   hev_.b = b;
0776   hev_.scale = scale;
0777   hev_.npart = npart;
0778   hev_.ncoll = ncoll;
0779   hev_.nhard = nhard;
0780   hev_.phi0 = phi0;
0781   hev_.vx = vx;
0782   hev_.vy = vy;
0783   hev_.vz = vz;
0784   hev_.vr = vr;
0785 
0786   if (doAnalysis_) {
0787     nt->Fill(nmix, np, src, sig);
0788     hydjetTree_->Fill();
0789   }
0790 
0791   //event counter
0792   if (doHistos_) {
0793     hNev->Fill(1., 1);
0794   }
0795 }
0796 // ------------ method called once each job just before starting event loop  ------------
0797 void HydjetAnalyzer::beginJob() {
0798   if (printLists_) {
0799     out_b.open(fBFileName.c_str());
0800     if (out_b.good() == false)
0801       throw cms::Exception("BadFile") << "Can\'t open file " << fBFileName;
0802     out_n.open(fNFileName.c_str());
0803     if (out_n.good() == false)
0804       throw cms::Exception("BadFile") << "Can\'t open file " << fNFileName;
0805     out_m.open(fMFileName.c_str());
0806     if (out_m.good() == false)
0807       throw cms::Exception("BadFile") << "Can\'t open file " << fMFileName;
0808   }
0809 
0810   if (doHistos_) {
0811     if (userHistos_) {
0812       //pt
0813       uhpt = new TH1D("uhpt", "uhpt", nintPt, ptBins);
0814       uhptj = new TH1D("uhptj", "uhptj", nintPt, ptBins);
0815       uhpth = new TH1D("uhpth", "uhpth", nintPt, ptBins);
0816 
0817       //pt_db
0818       uhpt_db = new TH1D("uhpt_db", "uhpt_db", 1000, 0.0000000000001, 100.);
0819       uhptj_db = new TH1D("uhptj_db", "uhptj_db", 1000, 0.0000000000001, 100.);
0820       uhpth_db = new TH1D("uhpth_db", "uhpth_db", 1000, 0.0000000000001, 100.);
0821 
0822       //eta
0823       uheta = new TH1D("uheta", "uheta", nintEta, etaBins);
0824       uhetaj = new TH1D("uhetaj", "uhetaj", nintEta, etaBins);
0825       uhetah = new TH1D("uhetah", "uhetah", nintEta, etaBins);
0826 
0827       //phi
0828       uhphi = new TH1D("uhphi", "uhphi", nintPhi, phiBins);
0829       uhphij = new TH1D("uhphij", "uhphij", nintPhi, phiBins);
0830       uhphih = new TH1D("uhphih", "uhphih", nintPhi, phiBins);
0831 
0832       const int NbinNpar = 5;
0833       const double BinsNpart[NbinNpar + 1] = {0., 29., 90., 202., 346., 400.};
0834 
0835       //ptNpart
0836       uhNpart = new TH1D("uhNpart", "uhNpart", NbinNpar, BinsNpart);
0837       uhNpartj = new TH1D("uhNpartj", "uhNpartj", NbinNpar, BinsNpart);
0838       uhNparth = new TH1D("uhNparth", "uhNparth", NbinNpar, BinsNpart);
0839 
0840       //ptNpart_db
0841       uhNpart_db = new TH1D("uhNpart_db", "uhNpart_db", 400, 0., 400.);
0842       uhNpartj_db = new TH1D("uhNpartj_db", "uhNpartj_db", 400, 0., 400.);
0843       uhNparth_db = new TH1D("uhNparth_db", "uhNparth_db", 400, 0., 400.);
0844 
0845       //ptNpart
0846       uhPtNpart = new TH1D("uhptNpart", "uhptNpart", NbinNpar, BinsNpart);
0847       uhPtNpartj = new TH1D("uhptNpartj", "uhptNpartj", NbinNpar, BinsNpart);
0848       uhPtNparth = new TH1D("uhptNparth", "uhptNparth", NbinNpar, BinsNpart);
0849 
0850       //ptNpart_db
0851       uhPtNpart_db = new TH1D("uhptNpart_db", "uhptNpart_db", 400, 0., 400.);
0852       uhPtNpartj_db = new TH1D("uhptNpartj_db", "uhptNpartj_db", 400, 0., 400.);
0853       uhPtNparth_db = new TH1D("uhptNparth_db", "uhptNparth_db", 400, 0., 400.);
0854 
0855       //v2Npart
0856       uhv2Npart = new TH1D("uhv2Npart", "uhv2Npart", NbinNpar, BinsNpart);
0857       uhv2Npartj = new TH1D("uhv2Npartj", "uhv2Npartj", NbinNpar, BinsNpart);
0858       uhv2Nparth = new TH1D("uhv2Nparth", "uhv2Nparth", NbinNpar, BinsNpart);
0859 
0860       //v2Npart_db
0861       uhv2Npart_db = new TH1D("uhv2Npart_db", "uhv2Npart_db", 400, 0., 400.);
0862       uhv2Npartj_db = new TH1D("uhv2Npartj_db", "uhv2Npartj_db", 400, 0., 400.);
0863       uhv2Nparth_db = new TH1D("uhv2Nparth_db", "uhv2Nparth_db", 400, 0., 400.);
0864 
0865       //v0pt
0866       uhv0pt = new TH1D("uhv0pt", "uhv0pt", nintV2pt, v2ptBins);
0867       uhv0ptj = new TH1D("uhv0ptj", "uhv0ptj", nintV2pt, v2ptBins);
0868       uhv0pth = new TH1D("uhv0pth", "uhv0pth", nintV2pt, v2ptBins);
0869 
0870       //v2pt
0871       uhv2pt = new TH1D("uhv2pt", "uhv2pt", nintV2pt, v2ptBins);
0872       uhv2ptj = new TH1D("uhv2ptj", "uhv2ptj", nintV2pt, v2ptBins);
0873       uhv2pth = new TH1D("uhv2pth", "uhv2pth", nintV2pt, v2ptBins);
0874 
0875       uhv3pt = new TH1D("uhv3pt", "uhv3pt", nintV2pt, v2ptBins);
0876       uhv4pt = new TH1D("uhv4pt", "uhv4pt", nintV2pt, v2ptBins);
0877       uhv5pt = new TH1D("uhv5pt", "uhv5pt", nintV2pt, v2ptBins);
0878       uhv6pt = new TH1D("uhv6pt", "uhv6pt", nintV2pt, v2ptBins);
0879 
0880       //v0pt
0881       uhv0pt_db = new TH1D("uhv0pt_db", "uhv0pt_db", 200, 0.0, 10.);
0882       uhv0ptj_db = new TH1D("uhv0ptj_db", "uhv0ptj_db", 200, 0.0, 10.);
0883       uhv0pth_db = new TH1D("uhv0pth_db", "uhv0pth_db", 200, 0.0, 10.);
0884 
0885       //v2pt_db
0886       uhv2pt_db = new TH1D("uhv2pt_db", "uhv2pt_db", 200, 0.0, 10.);
0887       uhv2ptj_db = new TH1D("uhv2ptj_db", "uhv2ptj_db", 200, 0.0, 10.);
0888       uhv2pth_db = new TH1D("uhv2pth_db", "uhv2pth_db", 200, 0.0, 10.);
0889 
0890       uhv3pt_db = new TH1D("uhv3pt_db", "uhv3pt_db", 200, 0.0, 10.);
0891       uhv4pt_db = new TH1D("uhv4pt_db", "uhv4pt_db", 200, 0.0, 10.);
0892       uhv5pt_db = new TH1D("uhv5pt_db", "uhv5pt_db", 200, 0.0, 10.);
0893       uhv6pt_db = new TH1D("uhv6pt_db", "uhv6pt_db", 200, 0.0, 10.);
0894 
0895       //v0eta
0896       uhv0eta = new TH1D("uhv0eta", "uhv0eta", nintV2eta, v2etaBins);
0897       uhv0etaj = new TH1D("uhv0etaj", "uhv0etaj", nintV2eta, v2etaBins);
0898       uhv0etah = new TH1D("uhv0etah", "uhv0etah", nintV2eta, v2etaBins);
0899 
0900       //v2eta
0901       uhv2eta = new TH1D("uhv2eta", "uhv2eta", nintV2eta, v2etaBins);
0902       uhv2etaj = new TH1D("uhv2etaj", "uhv2etaj", nintV2eta, v2etaBins);
0903       uhv2etah = new TH1D("uhv2etah", "uhv2etah", nintV2eta, v2etaBins);
0904 
0905       uhv3eta = new TH1D("uhv3eta", "uhv3eta", nintV2eta, v2etaBins);
0906       uhv4eta = new TH1D("uhv4eta", "uhv4eta", nintV2eta, v2etaBins);
0907       uhv5eta = new TH1D("uhv5eta", "uhv5eta", nintV2eta, v2etaBins);
0908       uhv6eta = new TH1D("uhv6eta", "uhv6eta", nintV2eta, v2etaBins);
0909 
0910       //v0eta_db
0911       uhv0eta_db = new TH1D("uhv0eta_db", "uhv0eta_db", 200, -5, 5.);
0912       uhv0etaj_db = new TH1D("uhv0etaj_db", "uhv0etaj_db", 200, -5, 5.);
0913       uhv0etah_db = new TH1D("uhv0etah_db", "uhv0etah_db", 200, -5, 5.);
0914 
0915       //v2eta_db
0916       uhv2eta_db = new TH1D("uhv2eta_db", "uhv2eta_db", 200, -5, 5.);
0917       uhv2etaj_db = new TH1D("uhv2etaj_db", "uhv2etaj_db", 200, -5, 5.);
0918       uhv2etah_db = new TH1D("uhv2etah_db", "uhv2etah_db", 200, -5, 5.);
0919 
0920       uhv3eta_db = new TH1D("uhv3eta_db", "uhv3eta_db", 200, -5, 5.);
0921       uhv4eta_db = new TH1D("uhv4eta_db", "uhv4eta_db", 200, -5, 5.);
0922       uhv5eta_db = new TH1D("uhv5eta_db", "uhv5eta_db", 200, -5, 5.);
0923       uhv6eta_db = new TH1D("uhv6eta_db", "uhv6eta_db", 200, -5, 5.);
0924     }
0925 
0926     dhpt = new TH1D("dhpt", "dhpt", 1000, 0.0000000000001, 100.);
0927     dhpt_ch = new TH1D("dhpt_ch", "dhpt_ch", 1000, 0.0000000000001, 100.);
0928 
0929     dhphi = new TH1D("dhphi", "dhphi", 1000, -3.14159265358979, 3.14159265358979);
0930     dhphi_ch = new TH1D("dhphi_ch", "dhphi_ch", 1000, -3.14159265358979, 3.14159265358979);
0931 
0932     dhv2pt = new TH1D("dhv2pt", "dhv2pt", 200, 0.0, 10.);
0933     dhv0pt = new TH1D("dhv0pt", "dhv0pt", 200, 0.0, 10.);
0934 
0935     dhv2eta = new TH1D("dhv2eta", "dhv2eta", 200, -5, 5.);
0936     dhv0eta = new TH1D("dhv0eta", "dhv0eta", 200, -5, 5.);
0937 
0938     dheta = new TH1D("dheta", "dheta", 1000, -10., 10.);
0939     dheta_ch = new TH1D("dheta_ch", "dheta_ch", 1000, -10., 10.);
0940 
0941     hNev = new TH1D("hNev", "hNev", 1, 0., 2.);
0942   }
0943 
0944   if (doAnalysis_) {
0945     usesResource(TFileService::kSharedResource);
0946     edm::Service<TFileService> f;
0947     nt = f->make<TNtuple>("nt", "Mixing Analysis", "mix:np:src:sig");
0948     hydjetTree_ = f->make<TTree>("hi", "Tree of Hydjet Events");
0949     hydjetTree_->Branch("event", &hev_.event, "event/I");
0950     hydjetTree_->Branch("b", &hev_.b, "b/F");
0951     hydjetTree_->Branch("npart", &hev_.npart, "npart/F");
0952     hydjetTree_->Branch("ncoll", &hev_.ncoll, "ncoll/F");
0953     hydjetTree_->Branch("nhard", &hev_.nhard, "nhard/F");
0954     hydjetTree_->Branch("phi0", &hev_.phi0, "phi0/F");
0955     hydjetTree_->Branch("scale", &hev_.scale, "scale/F");
0956     hydjetTree_->Branch("n", hev_.n, "n[3]/I");
0957     hydjetTree_->Branch("ptav", hev_.ptav, "ptav[3]/F");
0958     if (doParticles_) {
0959       hydjetTree_->Branch("mult", &hev_.mult, "mult/I");
0960       hydjetTree_->Branch("px", hev_.px, "px[mult]/F");
0961       hydjetTree_->Branch("py", hev_.py, "py[mult]/F");
0962       hydjetTree_->Branch("pz", hev_.pz, "pz[mult]/F");
0963       hydjetTree_->Branch("e", hev_.e, "e[mult]/F");
0964       hydjetTree_->Branch("pseudoRapidity", hev_.pseudoRapidity, "pseudoRapidity[mult]/F");
0965       hydjetTree_->Branch("pt", hev_.pt, "pt[mult]/F");
0966       hydjetTree_->Branch("eta", hev_.eta, "eta[mult]/F");
0967       hydjetTree_->Branch("phi", hev_.phi, "phi[mult]/F");
0968       hydjetTree_->Branch("pdg", hev_.pdg, "pdg[mult]/I");
0969       hydjetTree_->Branch("chg", hev_.chg, "chg[mult]/I");
0970 
0971       hydjetTree_->Branch("vx", &hev_.vx, "vx/F");
0972       hydjetTree_->Branch("vy", &hev_.vy, "vy/F");
0973       hydjetTree_->Branch("vz", &hev_.vz, "vz/F");
0974       hydjetTree_->Branch("vr", &hev_.vr, "vr/F");
0975     }
0976   }
0977 }
0978 // ------------ method called once each job just after ending the event loop  ------------
0979 void HydjetAnalyzer::endJob() {
0980   if (doHistos_) {
0981     dhpt->Write();
0982     dheta->Write();
0983     dhphi->Write();
0984     dhv0pt->Write();
0985     dhv2pt->Write();
0986     dhv0eta->Write();
0987     dhv2eta->Write();
0988 
0989     dhpt_ch->Write();
0990     dheta_ch->Write();
0991     hNev->Write();
0992     if (userHistos_) {
0993       uhpt->Write();
0994       uhpth->Write();
0995       uhptj->Write();
0996 
0997       uhpt_db->Write();
0998       uhpth_db->Write();
0999       uhptj_db->Write();
1000 
1001       uhNpart->Write();
1002       uhNparth->Write();
1003       uhNpartj->Write();
1004 
1005       uhNpart_db->Write();
1006       uhNparth_db->Write();
1007       uhNpartj_db->Write();
1008 
1009       uhPtNpart->Write();
1010       uhPtNparth->Write();
1011       uhPtNpartj->Write();
1012 
1013       uhPtNpart_db->Write();
1014       uhPtNparth_db->Write();
1015       uhPtNpartj_db->Write();
1016 
1017       uhv2Npart->Write();
1018       uhv2Nparth->Write();
1019       uhv2Npartj->Write();
1020 
1021       uhv2Npart_db->Write();
1022       uhv2Nparth_db->Write();
1023       uhv2Npartj_db->Write();
1024 
1025       uheta->Write();
1026       uhetah->Write();
1027       uhetaj->Write();
1028       uhphi->Write();
1029       uhphih->Write();
1030       uhphij->Write();
1031 
1032       uhv0eta->Write();
1033       uhv0etah->Write();
1034       uhv0etaj->Write();
1035 
1036       uhv0eta_db->Write();
1037       uhv0etah_db->Write();
1038       uhv0etaj_db->Write();
1039 
1040       uhv0pt->Write();
1041       uhv0pth->Write();
1042       uhv0ptj->Write();
1043 
1044       uhv0pt_db->Write();
1045       uhv0pth_db->Write();
1046       uhv0ptj_db->Write();
1047 
1048       uhv2eta->Write();
1049       uhv2etah->Write();
1050       uhv2etaj->Write();
1051 
1052       uhv2eta_db->Write();
1053       uhv2etah_db->Write();
1054       uhv2etaj_db->Write();
1055 
1056       uhv2pt->Write();
1057       uhv2pth->Write();
1058       uhv2ptj->Write();
1059 
1060       uhv2pt_db->Write();
1061       uhv2pth_db->Write();
1062       uhv2ptj_db->Write();
1063 
1064       uhv3eta->Write();
1065       uhv4eta->Write();
1066       uhv5eta->Write();
1067       uhv6eta->Write();
1068 
1069       uhv3eta_db->Write();
1070       uhv4eta_db->Write();
1071       uhv5eta_db->Write();
1072       uhv6eta_db->Write();
1073 
1074       uhv3pt->Write();
1075       uhv4pt->Write();
1076       uhv5pt->Write();
1077       uhv6pt->Write();
1078 
1079       uhv3pt_db->Write();
1080       uhv4pt_db->Write();
1081       uhv5pt_db->Write();
1082       uhv6pt_db->Write();
1083     }
1084   }
1085 }
1086 //define this as a plug-in
1087 DEFINE_FWK_MODULE(HydjetAnalyzer);