HydjetInterface/test/HydjetAnalyzer.cc

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