Hydjet2Interface/test/Hydjet2Analyzer.cc

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