Hydjet2Interface/test/Hydjet2Analyzer.cc

0001 /**
0002    \class Hydjet2Analyzer
0003    \brief HepMC events analyzer
0004    \version 2.2
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 #include "TH2.h"
0047
0048 #include "FWCore/Framework/interface/ConsumesCollector.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
0164   TH1D *dhphi;
0165
0166   TH1D *dhpdg;
0167
0168   TH1D *dhet_sum;
0169   TH1D *dhet_barrel_sum;
0170   TH1D *dhe_sum;
0171   TH1D *dhe_barrel_sum;
0172   TH1D *dheta;
0173   TH1D *dhpt;
0174
0175   TH1D *dhv2pt_cha;
0176   TH1D *dhv0pt_cha;
0177   TH1D *dhv2eta_cha;
0178   TH1D *dhv0eta_cha;
0179   TH1D *dhphi_cha;
0180   TH1D *dhet_cha_sum;
0181   TH1D *dhet_cha_barrel_sum;
0182   TH1D *dheta_cha;
0183   TH1D *dhpt_cha;
0184
0185   TH1D *dhv2pt_ch;
0186   TH1D *dhv0pt_ch;
0187   TH1D *dhv2eta_ch;
0188   TH1D *dhv0eta_ch;
0189   TH1D *dhphi_ch;
0190   TH1D *dhet_ch_sum;
0191   TH1D *dhet_ch_barrel_sum;
0192   TH1D *dhe_ch_sum;
0193   TH1D *dhe_ch_barrel_sum;
0194   TH1D *dheta_ch;
0195   TH1D *dhpt_ch;
0196
0197   TH1D *dhet_ph_sum;
0198   TH1D *dhet_ph_barrel_sum;
0199   TH1D *dhe_ph_sum;
0200   TH1D *dhe_ph_barrel_sum;
0201   TH1D *dheta_ph;
0202   TH1D *dhpt_ph;
0203
0204   TH1D *dhet_n_sum;
0205   TH1D *dhet_n_barrel_sum;
0206   TH1D *dhe_n_sum;
0207   TH1D *dhe_n_barrel_sum;
0208   TH1D *dheta_n;
0209   TH1D *dhpt_n;
0210
0211   TH1D *dhet_p_sum;
0212   TH1D *dhet_p_barrel_sum;
0213   TH1D *dheta_p;
0214   TH1D *dhpt_p;
0215
0216   TH1D *dhet_pi_sum;
0217   TH1D *dhet_pi_barrel_sum;
0218   TH1D *dheta_pi;
0219   TH1D *dhpt_pi;
0220
0221   TH1D *dhet_K_sum;
0222   TH1D *dhet_K_barrel_sum;
0223   TH1D *dheta_K;
0224   TH1D *dhpt_K;
0225
0226   TH2D *dhpdg_st;
0227
0228   TH1D *hNev;
0229
0230   //Users
0231   TH1D *uhpt;
0232   TH1D *uhpth;
0233   TH1D *uhptj;
0234
0235   TH1D *uhpt_db;
0236   TH1D *uhpth_db;
0237   TH1D *uhptj_db;
0238
0239   TH1D *uhNpart;
0240   TH1D *uhNparth;
0241   TH1D *uhNpartj;
0242
0243   TH1D *uhNpart_db;
0244   TH1D *uhNparth_db;
0245   TH1D *uhNpartj_db;
0246
0247   TH1D *uhPtNpart;
0248   TH1D *uhPtNparth;
0249   TH1D *uhPtNpartj;
0250
0251   TH1D *uhPtNpart_db;
0252   TH1D *uhPtNparth_db;
0253   TH1D *uhPtNpartj_db;
0254
0255   TH1D *uhv2Npart;
0256   TH1D *uhv2Nparth;
0257   TH1D *uhv2Npartj;
0258
0259   TH1D *uhv2Npart_db;
0260   TH1D *uhv2Nparth_db;
0261   TH1D *uhv2Npartj_db;
0262
0263   TH1D *uheta;
0264   TH1D *uhetah;
0265   TH1D *uhetaj;
0266
0267   TH1D *uhphi;
0268   TH1D *uhphih;
0269   TH1D *uhphij;
0270
0271   TH1D *uhv2pt;
0272   TH1D *uhv2pth;
0273   TH1D *uhv2ptj;
0274
0275   TH1D *uhv3pt;
0276   TH1D *uhv4pt;
0277   TH1D *uhv5pt;
0278   TH1D *uhv6pt;
0279
0280   TH1D *uhv0pt;
0281   TH1D *uhv0pth;
0282   TH1D *uhv0ptj;
0283
0284   TH1D *uhv2pt_db;
0285   TH1D *uhv2pth_db;
0286   TH1D *uhv2ptj_db;
0287
0288   TH1D *uhv3pt_db;
0289   TH1D *uhv4pt_db;
0290   TH1D *uhv5pt_db;
0291   TH1D *uhv6pt_db;
0292
0293   TH1D *uhv0pt_db;
0294   TH1D *uhv0pth_db;
0295   TH1D *uhv0ptj_db;
0296
0297   TH1D *uhv2eta;
0298   TH1D *uhv2etah;
0299   TH1D *uhv2etaj;
0300
0301   TH1D *uhv3eta;
0302   TH1D *uhv4eta;
0303   TH1D *uhv5eta;
0304   TH1D *uhv6eta;
0305
0306   TH1D *uhv0eta;
0307   TH1D *uhv0etah;
0308   TH1D *uhv0etaj;
0309
0310   TH1D *uhv2eta_db;
0311   TH1D *uhv2etah_db;
0312   TH1D *uhv2etaj_db;
0313
0314   TH1D *uhv3eta_db;
0315   TH1D *uhv4eta_db;
0316   TH1D *uhv5eta_db;
0317   TH1D *uhv6eta_db;
0318
0319   TH1D *uhv0eta_db;
0320   TH1D *uhv0etah_db;
0321   TH1D *uhv0etaj_db;
0322 };
0323
0324 Hydjet2Analyzer::Hydjet2Analyzer(const edm::ParameterSet &iConfig) {
0325   fBFileName = iConfig.getUntrackedParameter<std::string>("output_b", "b_values.txt");
0326   fNFileName = iConfig.getUntrackedParameter<std::string>("output_n", "n_values.txt");
0327   fMFileName = iConfig.getUntrackedParameter<std::string>("output_m", "m_values.txt");
0328   doAnalysis_ = iConfig.getUntrackedParameter<bool>("doAnalysis", false);
0329   useHepMCProduct_ = iConfig.getUntrackedParameter<bool>("useHepMCProduct", true);
0330   printLists_ = iConfig.getUntrackedParameter<bool>("printLists", false);
0331   doCF_ = iConfig.getUntrackedParameter<bool>("doMixed", false);
0332   doVertex_ = iConfig.getUntrackedParameter<bool>("doVertex", false);
0333   if (doVertex_) {
0334     simVerticesTag_ = iConfig.getParameter<edm::InputTag>("simVerticesTag");
0335   }
0336   etaMax_ = iConfig.getUntrackedParameter<double>("etaMax", 2.);
0337   ptMin_ = iConfig.getUntrackedParameter<double>("ptMin", 0);
0338   srcT_ = mayConsume<HepMCProduct>(
0339       iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("generator", "unsmeared")));
0340   srcTmix_ = consumes<CrossingFrame<edm::HepMCProduct>>(
0341       iConfig.getUntrackedParameter<edm::InputTag>("srcMix", edm::InputTag("mix", "generatorSmeared")));
0342
0343   genParticleSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("hiGenParticles"));
0344   genHIsrc_ = iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("heavyIon"));
0345
0346   if (useHepMCProduct_)
0347     pdtToken_ = esConsumes();
0348   doTestEvent_ = iConfig.getUntrackedParameter<bool>("doTestEvent", false);
0349   doParticles_ = iConfig.getUntrackedParameter<bool>("doParticles", false);
0350   doHistos_ = iConfig.getUntrackedParameter<bool>("doHistos", false);
0351   if (doHistos_) {
0352     userHistos_ = iConfig.getUntrackedParameter<bool>("userHistos", false);
0353     if (userHistos_) {
0354       uStatus_ = iConfig.getUntrackedParameter<int>("uStatus");
0355       uPDG_1 = iConfig.getUntrackedParameter<int>("uPDG_1");
0356       uPDG_2 = iConfig.getUntrackedParameter<int>("uPDG_2", uPDG_1);
0357       uPDG_3 = iConfig.getUntrackedParameter<int>("uPDG_3", uPDG_1);
0358       upTetaCut_ = iConfig.getUntrackedParameter<double>("uPTetaCut", 0.8);
0359       downTetaCut_ = iConfig.getUntrackedParameter<double>("dPTetaCut", -1.);
0360       uPtBins_ = iConfig.getUntrackedParameter<vector<double>>("PtBins");
0361       uEtaBins_ = iConfig.getUntrackedParameter<vector<double>>("EtaBins");
0362       uPhiBins_ = iConfig.getUntrackedParameter<vector<double>>("PhiBins");
0363       uV2ptBins_ = iConfig.getUntrackedParameter<vector<double>>("v2PtBins");
0364       uV2etaBins_ = iConfig.getUntrackedParameter<vector<double>>("v2EtaBins");
0365
0366       //Pt
0367       int PtSize = uPtBins_.size();
0368       if (PtSize > 1) {
0369         ptBins = new float[PtSize];
0370         nintPt = PtSize - 1;
0371         for (int k = 0; k < PtSize; k++) {
0372           ptBins[k] = uPtBins_[k];
0373         }
0374       } else {
0375         nintPt = iConfig.getUntrackedParameter<int>("nintPt");
0376         maxPt = iConfig.getUntrackedParameter<double>("maxPt");
0377         minPt = iConfig.getUntrackedParameter<double>("minPt");
0378         ptBins = new float[nintPt + 1];
0379         for (int k = 0; k < nintPt + 1; k++) {
0380           ptBins[k] = minPt + k * ((maxPt - minPt) / nintPt);
0381         }
0382       }
0383
0384       //Eta
0385       int EtaSize = uEtaBins_.size();
0386       if (EtaSize > 1) {
0387         etaBins = new float[EtaSize];
0388         nintEta = EtaSize - 1;
0389         for (int k = 0; k < EtaSize; k++) {
0390           etaBins[k] = uEtaBins_[k];
0391         }
0392       } else {
0393         nintEta = iConfig.getUntrackedParameter<int>("nintEta");
0394         maxEta = iConfig.getUntrackedParameter<double>("maxEta");
0395         minEta = iConfig.getUntrackedParameter<double>("minEta");
0396         etaBins = new float[nintEta + 1];
0397         for (int k = 0; k < nintEta + 1; k++) {
0398           etaBins[k] = minEta + k * ((maxEta - minEta) / nintEta);
0399         }
0400       }
0401
0402       //Phi
0403       int PhiSize = uPhiBins_.size();
0404       if (PhiSize > 1) {
0405         phiBins = new float[PhiSize];
0406         nintPhi = PhiSize - 1;
0407         for (int k = 0; k < PhiSize; k++) {
0408           phiBins[k] = uPhiBins_[k];
0409         }
0410       } else {
0411         nintPhi = iConfig.getUntrackedParameter<int>("nintPhi");
0412         maxPhi = iConfig.getUntrackedParameter<double>("maxPhi");
0413         minPhi = iConfig.getUntrackedParameter<double>("minPhi");
0414         phiBins = new float[nintPhi + 1];
0415         for (int k = 0; k < nintPhi + 1; k++) {
0416           phiBins[k] = minPhi + k * ((maxPhi - minPhi) / nintPhi);
0417         }
0418       }
0419
0420       //v2Pt
0421       int v2PtSize = uV2ptBins_.size();
0422       if (v2PtSize > 1) {
0423         v2ptBins = new float[v2PtSize];
0424         nintV2pt = v2PtSize - 1;
0425         for (int k = 0; k < v2PtSize; k++) {
0426           v2ptBins[k] = uV2ptBins_[k];
0427         }
0428       } else {
0429         nintV2pt = iConfig.getUntrackedParameter<int>("nintV2pt");
0430         maxV2pt = iConfig.getUntrackedParameter<double>("maxV2pt");
0431         minV2pt = iConfig.getUntrackedParameter<double>("minV2pt");
0432         v2ptBins = new float[nintV2pt + 1];
0433         for (int k = 0; k < nintV2pt + 1; k++) {
0434           v2ptBins[k] = minV2pt + k * ((maxV2pt - minV2pt) / nintV2pt);
0435         }
0436       }
0437
0438       //v2Eta
0439       int v2EtaSize = uV2etaBins_.size();
0440       if (v2EtaSize > 1) {
0441         v2etaBins = new float[v2EtaSize];
0442         nintV2eta = v2EtaSize - 1;
0443         for (int k = 0; k < v2EtaSize; k++) {
0444           v2etaBins[k] = uV2etaBins_[k];
0445         }
0446       } else {
0447         nintV2eta = iConfig.getUntrackedParameter<int>("nintV2eta");
0448         maxV2eta = iConfig.getUntrackedParameter<double>("maxV2eta");
0449         minV2eta = iConfig.getUntrackedParameter<double>("minV2eta");
0450         v2etaBins = new float[nintV2eta + 1];
0451         for (int k = 0; k < nintV2eta + 1; k++) {
0452           v2etaBins[k] = minV2eta + k * ((maxV2eta - minV2eta) / nintV2eta);
0453         }
0454       }
0455     }  //user histo
0456   }  //do histo
0457 }
0458 Hydjet2Analyzer::~Hydjet2Analyzer() {}
0459
0460 // ------------ method called to for each event  ------------
0461 void Hydjet2Analyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0462   using namespace edm;
0463   using namespace HepMC;
0464   hev_.event = iEvent.id().event();
0465   for (int ieta = 0; ieta < ETABINS; ++ieta) {
0466     hev_.n[ieta] = 0;
0467     hev_.ptav[ieta] = 0;
0468   }
0469   hev_.mult = 0;
0470   double phi0 = 0.;
0471   double phi3 = 0.;
0472   double b = -1.;
0473   double v2, v3, v4, v5, v6;
0474   double scale = -1.;
0475   int npart = -1;
0476   int ncoll = -1;
0477   int nhard = -1;
0478   double vx = -99.;
0479   double vy = -99.;
0480   double vz = -99.;
0481   double vr = -99.;
0482   const GenEvent *evt;
0483   int nmix = -1;
0484   int np = 0;
0485   int sig = -1;
0486   int src = -1;
0487   if (useHepMCProduct_) {
0488     edm::ESHandle<ParticleDataTable> pdt = iSetup.getHandle(pdtToken_);
0489
0490     if (doCF_) {
0491       Handle<CrossingFrame<HepMCProduct>> cf;
0492       iEvent.getByToken(srcTmix_, cf);
0493       MixCollection<HepMCProduct> mix(cf.product());
0494       nmix = mix.size();
0495       //cout << "Mix Collection Size: " << mix.size() <<", pileup size: " <<mix.sizePileup() << ", signal: "<<mix.sizeSignal()<< endl;
0496       MixCollection<HepMCProduct>::iterator mbegin = mix.begin();
0497       MixCollection<HepMCProduct>::iterator mend = mix.end();
0498       for (MixCollection<HepMCProduct>::iterator mixit = mbegin; mixit != mend; ++mixit) {
0499         const GenEvent *subevt = (*mixit).GetEvent();
0500         int all = subevt->particles_size();
0501         //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;
0502         np += all;
0503         HepMC::GenEvent::particle_const_iterator begin = subevt->particles_begin();
0504         HepMC::GenEvent::particle_const_iterator end = subevt->particles_end();
0505         for (HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it) {
0506           if ((*it)->status() == 1) {
0507             int pdg_id = (*it)->pdg_id();
0508             float eta = (*it)->momentum().eta();
0509             float phi = (*it)->momentum().phi();
0510             float pt = (*it)->momentum().perp();
0511             const ParticleData *part = pdt->particle(pdg_id);
0512             int charge = static_cast<int>(part->charge());
0513             hev_.pt[hev_.mult] = pt;
0514             hev_.eta[hev_.mult] = eta;
0515             hev_.phi[hev_.mult] = phi;
0516             hev_.pdg[hev_.mult] = pdg_id;
0517             hev_.chg[hev_.mult] = charge;
0518
0519             //cout << "Mix Particles: pt= " << pt<<" eta="<<eta<<" phi="<<phi<<" pdg="<< pdg_id<<" charge="<<charge << endl;
0520             eta = fabs(eta);
0521             int etabin = 0;
0522             if (eta > 0.5)
0523               etabin = 1;
0524             if (eta > 1.)
0525               etabin = 2;
0526             if (eta < 2.) {
0527               hev_.ptav[etabin] += pt;
0528               ++(hev_.n[etabin]);
0529             }
0530             ++(hev_.mult);
0531           }
0532         }
0533       }
0534     } else {  //not mixing
0535       Handle<HepMCProduct> mc;
0536       iEvent.getByToken(srcT_, mc);
0537       evt = mc->GetEvent();
0538       scale = evt->event_scale();
0539       const HeavyIon *hi = evt->heavy_ion();
0540
0541       if (hi) {
0542         b = hi->impact_parameter();
0543         npart = hi->Npart_proj() + hi->Npart_targ();
0544         ncoll = hi->Ncoll();
0545         nhard = hi->Ncoll_hard();
0546         phi0 = hi->event_plane_angle();
0547         phi3 =
0548             hi->eccentricity();  // 0.;  // No HepMC entry for Psi3 exist, but in private code it's possible to use hi->eccentricity();
0549         if (printLists_) {
0550           out_b << b << endl;
0551           out_n << npart << endl;
0552         }
0553       }
0554
0555       src = evt->particles_size();
0556       if (doTestEvent_) {
0557         std::cout << "Event size: " << src << std::endl;
0558         mc->GetEvent()->print();
0559       }
0560
0561       float et_sum = 0., et_barrel_sum = 0., e_sum = 0., e_barrel_sum = 0., et_cha_sum = 0., et_cha_barrel_sum = 0.,
0562             et_ch_sum = 0., et_ch_barrel_sum = 0., e_ch_sum = 0., e_ch_barrel_sum = 0., et_ph_sum = 0.,
0563             et_ph_barrel_sum = 0., e_ph_sum = 0., e_ph_barrel_sum = 0., et_n_sum = 0., et_n_barrel_sum = 0.,
0564             et_p_sum = 0., et_p_barrel_sum = 0., et_pi_sum = 0., et_pi_barrel_sum = 0., et_K_sum = 0.,
0565             et_K_barrel_sum = 0., e_n_sum = 0., e_n_barrel_sum = 0.;
0566
0567       HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
0568       HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
0569       for (HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it) {
0570         if (((*it)->status() >= -1) && ((*it)->status() < 31)) {
0571           const ParticleData *part;
0572           int st = (*it)->status();
0573           int pdg_id = (*it)->pdg_id();
0574           float eta = (*it)->momentum().eta();
0575           float phi = (*it)->momentum().phi();
0576           float pt = (*it)->momentum().perp();
0577
0578           float px = (*it)->momentum().px();
0579           float py = (*it)->momentum().py();
0580           float pz = (*it)->momentum().pz();
0581           float e = (*it)->momentum().e();
0582           float pseudoRapidity = (*it)->momentum().pseudoRapidity();
0583           int charge = -1;
0584           float mass = -1.;
0585           if ((pdt->particle(pdg_id))) {
0586             part = pdt->particle(pdg_id);
0587             charge = static_cast<int>(part->charge());
0588             mass = (part->mass());
0589           }
0590           if (mass < 0) {
0591             if ((abs(pdg_id) == 130) || (abs(pdg_id) == 310)) {
0592               mass = 0.497611;
0593             } else {
0594               //cout<<"Error! Mass for "<< pdg_id <<" not found in PDT!!!"<<endl;
0595               //return;
0596             }
0597           }
0598
0599           float et = sqrt((pt * pt) + (mass * mass));
0600
0601           //if(std::abs(eta)>5) continue;
0602
0603           dhpdg_st->Fill(pdg_id, st);
0604
0605           hev_.px[hev_.mult] = px;
0606           hev_.py[hev_.mult] = py;
0607           hev_.pz[hev_.mult] = pz;
0608           hev_.e[hev_.mult] = e;
0609           hev_.pseudoRapidity[hev_.mult] = pseudoRapidity;
0610           hev_.pt[hev_.mult] = pt;
0611           hev_.eta[hev_.mult] = eta;
0612           hev_.phi[hev_.mult] = phi;
0613           hev_.pdg[hev_.mult] = pdg_id;
0614           hev_.chg[hev_.mult] = charge;
0615
0616           phi = phi - phi0;
0617           ///
0618           double phiTrue;
0619           if (phi > pi) {
0620             phiTrue = phi - (2 * pi);
0621           } else if (phi < (-1 * pi)) {
0622             phiTrue = phi + (2 * pi);
0623           } else {
0624             phiTrue = phi;
0625           }
0626           ///
0627           v2 = std::cos(2 * (phiTrue));
0628           v3 = std::cos(3 * (phiTrue - phi3));
0629           v4 = std::cos(4 * (phiTrue));
0630           v5 = std::cos(5 * (phiTrue - phi3));
0631           v6 = std::cos(6 * (phiTrue));
0632
0633           if (doHistos_) {
0634             //common histos
0635             if ((*it)->status() == 1) {
0636               dhpdg->Fill(pdg_id);
0637               dheta->Fill(eta);
0638
0639               e_sum += e;
0640               et_sum += et;
0641               if (std::abs(eta) < 1.5) {
0642                 e_barrel_sum += e;
0643                 et_barrel_sum += et;
0644               }
0645
0646               if (std::abs(pdg_id) == 22) {  //ph
0647                 e_ph_sum += e;
0648                 et_ph_sum += et;
0649                 if (std::abs(eta) < 1.5) {
0650                   e_ph_barrel_sum += e;
0651                   et_ph_barrel_sum += et;
0652                 }
0653                 if (std::abs(eta) < 0.8) {  //ALICE
0654                   dhpt_ph->Fill(pt);
0655                 }
0656                 dheta_ph->Fill(eta);
0657               }
0658
0659               if (std::abs(pdg_id) == 2112) {  //n
0660                 e_n_sum += e;
0661                 et_n_sum += et;
0662                 if (std::abs(eta) < 1.5) {
0663                   e_n_barrel_sum += e;
0664                   et_n_barrel_sum += et;
0665                 }
0666                 if (std::abs(eta) < 0.8) {
0667                   dhpt_n->Fill(pt);
0668                 }
0669                 dheta_n->Fill(eta);
0670               }
0671
0672               if (std::abs(pdg_id) == 2212) {  //p
0673                 et_p_sum += et;
0674                 if (std::abs(eta) < 1.5) {
0675                   et_p_barrel_sum += et;
0676                 }
0677                 if (std::abs(eta) < 0.8) {
0678                   dhpt_p->Fill(pt);
0679                 }
0680                 dheta_p->Fill(eta);
0681               }
0682
0683               if (std::abs(pdg_id) == 211) {  //pi
0684                 et_pi_sum += et;
0685                 if (std::abs(eta) < 1.5) {
0686                   et_pi_barrel_sum += et;
0687                 }
0688                 if (std::abs(eta) < 0.8) {
0689                   dhpt_pi->Fill(pt);
0690                 }
0691                 dheta_pi->Fill(eta);
0692               }
0693
0694               if (std::abs(pdg_id) == 321) {  //K
0695                 et_K_sum += et;
0696                 if (std::abs(eta) < 1.5) {
0697                   et_K_barrel_sum += et;
0698                 }
0699                 if (std::abs(eta) < 0.8) {
0700                   dhpt_K->Fill(pt);
0701                 }
0702                 dheta_K->Fill(eta);
0703               }
0704
0705               if (std::abs(eta) < 0.8) {
0706                 dhpt->Fill(pt);
0707                 dhphi->Fill(phiTrue);
0708               }
0709
0710               if (charge == 1) {
0711                 et_cha_sum += et;
0712                 if (std::abs(eta) < 1.5) {
0713                   et_cha_barrel_sum += et;
0714                 }
0715
0716                 if (std::abs(eta) < 1.) {  //CMS
0717                   dhv0pt_cha->Fill(pt, 1.);
0718                   dhv2pt_cha->Fill(pt, v2);
0719                 }
0720
0721                 if (std::abs(eta) < 0.8) {  //ALICE
0722                   dhpt_cha->Fill(pt);
0723                   dhphi_cha->Fill(phiTrue);
0724                 }
0725                 dhv0eta_cha->Fill(eta, 1.);
0726                 dhv2eta_cha->Fill(eta, v2);
0727                 dheta_cha->Fill(eta);
0728               }
0729
0730               if (std::abs(pdg_id) == 211 || std::abs(pdg_id) == 321 || std::abs(pdg_id) == 2212) {  //ch
0731                 et_ch_sum += et;
0732                 e_ch_sum += e;
0733
0734                 if (std::abs(eta) < 1.5) {
0735                   e_ch_barrel_sum += e;
0736                   et_ch_barrel_sum += et;
0737                 }
0738
0739                 if (std::abs(eta) < 0.8) {
0740                   dhv0pt_ch->Fill(pt, 1.);
0741                   dhv2pt_ch->Fill(pt, v2);
0742                   dhpt_ch->Fill(pt);
0743                   dhphi_ch->Fill(phiTrue);
0744                 }
0745
0746                 dhv0eta_ch->Fill(eta, 1.);
0747                 dhv2eta_ch->Fill(eta, v2);
0748                 dheta_ch->Fill(eta);
0749               }  //ch
0750             }  //status 1
0751
0752             //user histos
0753             if (userHistos_ && ((uStatus_ == 3) || (((*it)->status() < 10) && (uStatus_ == 1)) ||
0754                                 (((*it)->status() > 10) && (uStatus_ == 2)))) {  //user status
0755
0756               //set1
0757               if (std::abs(pdg_id) == uPDG_1 || std::abs(pdg_id) == uPDG_2 || std::abs(pdg_id) == uPDG_3) {  //uPDG
0758                 if ((uStatus_ == 3) && ((*it)->status() < 10))
0759                   cout << "ustatus=3, but stab. part. found!!!" << endl;
0760
0761                 if (std::abs(eta) > downTetaCut_ && std::abs(eta) < upTetaCut_) {  //eta cut
0762
0763                   uhv0pt->Fill(pt, 1.);
0764                   uhv2pt->Fill(pt, v2);
0765                   uhv3pt->Fill(pt, v3);
0766                   uhv4pt->Fill(pt, v4);
0767                   uhv5pt->Fill(pt, v5);
0768                   uhv6pt->Fill(pt, v6);
0769
0770                   uhv0pt_db->Fill(pt, 1.);
0771                   uhv2pt_db->Fill(pt, v2);
0772                   uhv3pt_db->Fill(pt, v3);
0773                   uhv4pt_db->Fill(pt, v4);
0774                   uhv5pt_db->Fill(pt, v5);
0775                   uhv6pt_db->Fill(pt, v6);
0776
0777                   if (pt >= 1.5 && pt < 10.) {
0778                     uhv2Npart->Fill(npart, v2);
0779                     uhv2Npart_db->Fill(npart, v2);
0780
0781                     uhPtNpart->Fill(npart, pt);
0782                     uhPtNpart_db->Fill(npart, pt);
0783
0784                     uhNpart->Fill(npart, 1.);
0785                     uhNpart_db->Fill(npart, 1.);
0786                   }
0787
0788                   uhpt->Fill(pt);
0789                   uhpt_db->Fill(pt);
0790                   uhphi->Fill(phiTrue);
0791
0792                   if (((*it)->status() == 16) || ((*it)->status() == 6)) {  //hydro
0793                     uhv0pth->Fill(pt, 1.);
0794                     uhv2pth->Fill(pt, v2);
0795
0796                     uhv0pth_db->Fill(pt, 1.);
0797                     uhv2pth_db->Fill(pt, v2);
0798
0799                     if (pt >= 1.5 && pt < 10.) {
0800                       uhv2Nparth->Fill(npart, v2);
0801                       uhv2Nparth_db->Fill(npart, v2);
0802                     }
0803
0804                     uhPtNparth->Fill(npart, pt);
0805                     uhPtNparth_db->Fill(npart, pt);
0806
0807                     uhpth->Fill(pt);
0808                     uhpth_db->Fill(pt);
0809                     uhphih->Fill(phiTrue);
0810                   }
0811
0812                   if (((*it)->status() == 17) || ((*it)->status() == 7)) {  //jet
0813                     uhv0ptj->Fill(pt, 1.);
0814                     uhv2ptj->Fill(pt, v2);
0815
0816                     uhv0ptj_db->Fill(pt, 1.);
0817                     uhv2ptj_db->Fill(pt, v2);
0818
0819                     if (pt >= 1.5 && pt < 10.) {
0820                       uhv2Npartj->Fill(npart, v2);
0821                       uhv2Npartj_db->Fill(npart, v2);
0822                     }
0823
0824                     uhPtNpartj->Fill(npart, pt);
0825                     uhPtNpartj_db->Fill(npart, pt);
0826
0827                     uhptj->Fill(pt);
0828                     uhptj_db->Fill(pt);
0829                     uhphij->Fill(phiTrue);
0830                   }
0831                 }  //eta cut
0832
0833                 uheta->Fill(eta);
0834
0835                 uhv0eta->Fill(eta, 1.);
0836                 uhv2eta->Fill(eta, v2);
0837                 uhv3eta->Fill(eta, v3);
0838                 uhv4eta->Fill(eta, v4);
0839                 uhv5eta->Fill(eta, v5);
0840                 uhv6eta->Fill(eta, v6);
0841
0842                 uhv0eta_db->Fill(eta, 1.);
0843                 uhv2eta_db->Fill(eta, v2);
0844                 uhv3eta_db->Fill(eta, v3);
0845                 uhv4eta_db->Fill(eta, v4);
0846                 uhv5eta_db->Fill(eta, v5);
0847                 uhv6eta_db->Fill(eta, v6);
0848
0849                 if (((*it)->status() == 16) || ((*it)->status() == 6)) {  //hydro
0850                   uhv2etah->Fill(eta, v2);
0851                   uhv0etah->Fill(eta, 1.);
0852
0853                   uhv2etah_db->Fill(eta, v2);
0854                   uhv0etah_db->Fill(eta, 1.);
0855
0856                   uhetah->Fill(eta);
0857                 }
0858                 if (((*it)->status() == 17) || ((*it)->status() == 7)) {  //jet
0859                   uhv2etaj->Fill(eta, v2);
0860                   uhv0etaj->Fill(eta, 1.);
0861
0862                   uhv2etaj_db->Fill(eta, v2);
0863                   uhv0etaj_db->Fill(eta, 1.);
0864
0865                   uhetaj->Fill(eta);
0866                 }
0867               }  //uPDG
0868
0869             }  //user status
0870
0871           }  //doHistos_
0872
0873           eta = fabs(eta);
0874           int etabin = 0;
0875           if (eta > 0.5)
0876             etabin = 1;
0877           if (eta > 1.)
0878             etabin = 2;
0879           if (eta < 2.) {
0880             hev_.ptav[etabin] += pt;
0881             ++(hev_.n[etabin]);
0882           }
0883           ++(hev_.mult);
0884         }
0885       }  //particle iterator
0886       dhet_sum->Fill(et_sum);
0887       dhet_barrel_sum->Fill(et_barrel_sum);
0888       dhe_sum->Fill(e_sum);
0889       dhe_barrel_sum->Fill(e_barrel_sum);
0890
0891       dhet_cha_sum->Fill(et_cha_sum);
0892       dhet_cha_barrel_sum->Fill(et_cha_barrel_sum);
0893
0894       dhet_ch_sum->Fill(et_ch_sum);
0895       dhet_ch_barrel_sum->Fill(et_ch_barrel_sum);
0896       dhe_ch_sum->Fill(e_ch_sum);
0897       dhe_ch_barrel_sum->Fill(e_ch_barrel_sum);
0898
0899       dhet_ph_sum->Fill(et_ph_sum);
0900       dhet_ph_barrel_sum->Fill(et_ph_barrel_sum);
0901       dhe_ph_sum->Fill(e_ph_sum);
0902       dhe_ph_barrel_sum->Fill(e_ph_barrel_sum);
0903
0904       dhet_n_sum->Fill(et_n_sum);
0905       dhet_n_barrel_sum->Fill(et_n_barrel_sum);
0906       dhe_n_sum->Fill(e_n_sum);
0907       dhe_n_barrel_sum->Fill(e_n_barrel_sum);
0908
0909       dhet_p_sum->Fill(et_p_sum);
0910       dhet_p_barrel_sum->Fill(et_p_barrel_sum);
0911       dhet_pi_sum->Fill(et_pi_sum);
0912       dhet_pi_barrel_sum->Fill(et_pi_barrel_sum);
0913       dhet_K_sum->Fill(et_K_sum);
0914       dhet_K_barrel_sum->Fill(et_K_barrel_sum);
0915
0916     }  //not mixing
0917   } else {  // not HepMC
0918     edm::Handle<reco::GenParticleCollection> parts;
0919     iEvent.getByLabel(genParticleSrc_, parts);
0920     for (unsigned int i = 0; i < parts->size(); ++i) {
0921       const reco::GenParticle &p = (*parts)[i];
0922       hev_.pt[hev_.mult] = p.pt();
0923       hev_.eta[hev_.mult] = p.eta();
0924       hev_.phi[hev_.mult] = p.phi();
0925       hev_.pdg[hev_.mult] = p.pdgId();
0926       hev_.chg[hev_.mult] = p.charge();
0927       double eta = fabs(p.eta());
0928
0929       int etabin = 0;
0930       if (eta > 0.5)
0931         etabin = 1;
0932       if (eta > 1.)
0933         etabin = 2;
0934       if (eta < 2.) {
0935         hev_.ptav[etabin] += p.pt();
0936         ++(hev_.n[etabin]);
0937       }
0938       ++(hev_.mult);
0939     }
0940     if (doHI_) {
0941       edm::Handle<GenHIEvent> higen;
0942       iEvent.getByLabel(genHIsrc_, higen);
0943     }
0944   }
0945
0946   if (doVertex_) {
0947     edm::Handle<edm::SimVertexContainer> simVertices;
0948     iEvent.getByLabel<edm::SimVertexContainer>(simVerticesTag_, simVertices);
0949
0950     if (!simVertices.isValid())
0951       throw cms::Exception("FatalError") << "No vertices found\n";
0952     int inum = 0;
0953
0954     edm::SimVertexContainer::const_iterator it = simVertices->begin();
0955     SimVertex vertex = (*it);
0956     cout << " Vertex position " << inum << " " << vertex.position().rho() << " " << vertex.position().z() << endl;
0957     vx = vertex.position().x();
0958     vy = vertex.position().y();
0959     vz = vertex.position().z();
0960     vr = vertex.position().rho();
0961   }
0962
0963   for (int i = 0; i < 3; ++i) {
0964     hev_.ptav[i] = hev_.ptav[i] / hev_.n[i];
0965   }
0966
0967   hev_.b = b;
0968   hev_.scale = scale;
0969   hev_.npart = npart;
0970   hev_.ncoll = ncoll;
0971   hev_.nhard = nhard;
0972   hev_.phi0 = phi0;
0973   hev_.vx = vx;
0974   hev_.vy = vy;
0975   hev_.vz = vz;
0976   hev_.vr = vr;
0977
0978   if (doAnalysis_) {
0979     nt->Fill(nmix, np, src, sig);
0980     hydjetTree_->Fill();
0981   }
0982
0983   //event counter
0984   if (doHistos_) {
0985     hNev->Fill(1., 1);
0986   }
0987 }
0988 // ------------ method called once each job just before starting event loop  ------------
0989 void Hydjet2Analyzer::beginJob() {
0990   if (printLists_) {
0991     out_b.open(fBFileName.c_str());
0992     if (out_b.good() == false)
0993       throw cms::Exception("BadFile") << "Can\'t open file " << fBFileName;
0994     out_n.open(fNFileName.c_str());
0995     if (out_n.good() == false)
0996       throw cms::Exception("BadFile") << "Can\'t open file " << fNFileName;
0997     out_m.open(fMFileName.c_str());
0998     if (out_m.good() == false)
0999       throw cms::Exception("BadFile") << "Can\'t open file " << fMFileName;
1000   }
1001
1002   if (doHistos_) {
1003     if (userHistos_) {
1004       //pt
1005       uhpt = new TH1D("uhpt", "uhpt", nintPt, ptBins);
1006       uhptj = new TH1D("uhptj", "uhptj", nintPt, ptBins);
1007       uhpth = new TH1D("uhpth", "uhpth", nintPt, ptBins);
1008
1009       //pt_db
1010       uhpt_db = new TH1D("uhpt_db", "uhpt_db", 1000, 0.0000000000001, 100.);
1011       uhptj_db = new TH1D("uhptj_db", "uhptj_db", 1000, 0.0000000000001, 100.);
1012       uhpth_db = new TH1D("uhpth_db", "uhpth_db", 1000, 0.0000000000001, 100.);
1013
1014       //eta
1015       uheta = new TH1D("uheta", "uheta", nintEta, etaBins);
1016       uhetaj = new TH1D("uhetaj", "uhetaj", nintEta, etaBins);
1017       uhetah = new TH1D("uhetah", "uhetah", nintEta, etaBins);
1018
1019       //phi
1020       uhphi = new TH1D("uhphi", "uhphi", nintPhi, phiBins);
1021       uhphij = new TH1D("uhphij", "uhphij", nintPhi, phiBins);
1022       uhphih = new TH1D("uhphih", "uhphih", nintPhi, phiBins);
1023
1024       const int NbinNpar = 5;
1025       const double BinsNpart[NbinNpar + 1] = {0., 29., 90., 202., 346., 400.};
1026
1027       //ptNpart
1028       uhNpart = new TH1D("uhNpart", "uhNpart", NbinNpar, BinsNpart);
1029       uhNpartj = new TH1D("uhNpartj", "uhNpartj", NbinNpar, BinsNpart);
1030       uhNparth = new TH1D("uhNparth", "uhNparth", NbinNpar, BinsNpart);
1031
1032       //ptNpart_db
1033       uhNpart_db = new TH1D("uhNpart_db", "uhNpart_db", 400, 0., 400.);
1034       uhNpartj_db = new TH1D("uhNpartj_db", "uhNpartj_db", 400, 0., 400.);
1035       uhNparth_db = new TH1D("uhNparth_db", "uhNparth_db", 400, 0., 400.);
1036
1037       //ptNpart
1038       uhPtNpart = new TH1D("uhptNpart", "uhptNpart", NbinNpar, BinsNpart);
1039       uhPtNpartj = new TH1D("uhptNpartj", "uhptNpartj", NbinNpar, BinsNpart);
1040       uhPtNparth = new TH1D("uhptNparth", "uhptNparth", NbinNpar, BinsNpart);
1041
1042       //ptNpart_db
1043       uhPtNpart_db = new TH1D("uhptNpart_db", "uhptNpart_db", 400, 0., 400.);
1044       uhPtNpartj_db = new TH1D("uhptNpartj_db", "uhptNpartj_db", 400, 0., 400.);
1045       uhPtNparth_db = new TH1D("uhptNparth_db", "uhptNparth_db", 400, 0., 400.);
1046
1047       //v2Npart
1048       uhv2Npart = new TH1D("uhv2Npart", "uhv2Npart", NbinNpar, BinsNpart);
1049       uhv2Npartj = new TH1D("uhv2Npartj", "uhv2Npartj", NbinNpar, BinsNpart);
1050       uhv2Nparth = new TH1D("uhv2Nparth", "uhv2Nparth", NbinNpar, BinsNpart);
1051
1052       //v2Npart_db
1053       uhv2Npart_db = new TH1D("uhv2Npart_db", "uhv2Npart_db", 400, 0., 400.);
1054       uhv2Npartj_db = new TH1D("uhv2Npartj_db", "uhv2Npartj_db", 400, 0., 400.);
1055       uhv2Nparth_db = new TH1D("uhv2Nparth_db", "uhv2Nparth_db", 400, 0., 400.);
1056
1057       //v0pt
1058       uhv0pt = new TH1D("uhv0pt", "uhv0pt", nintV2pt, v2ptBins);
1059       uhv0ptj = new TH1D("uhv0ptj", "uhv0ptj", nintV2pt, v2ptBins);
1060       uhv0pth = new TH1D("uhv0pth", "uhv0pth", nintV2pt, v2ptBins);
1061
1062       //v2pt
1063       uhv2pt = new TH1D("uhv2pt", "uhv2pt", nintV2pt, v2ptBins);
1064       uhv2ptj = new TH1D("uhv2ptj", "uhv2ptj", nintV2pt, v2ptBins);
1065       uhv2pth = new TH1D("uhv2pth", "uhv2pth", nintV2pt, v2ptBins);
1066
1067       uhv3pt = new TH1D("uhv3pt", "uhv3pt", nintV2pt, v2ptBins);
1068       uhv4pt = new TH1D("uhv4pt", "uhv4pt", nintV2pt, v2ptBins);
1069       uhv5pt = new TH1D("uhv5pt", "uhv5pt", nintV2pt, v2ptBins);
1070       uhv6pt = new TH1D("uhv6pt", "uhv6pt", nintV2pt, v2ptBins);
1071
1072       //v0pt
1073       uhv0pt_db = new TH1D("uhv0pt_db", "uhv0pt_db", 200, 0.0, 10.);
1074       uhv0ptj_db = new TH1D("uhv0ptj_db", "uhv0ptj_db", 200, 0.0, 10.);
1075       uhv0pth_db = new TH1D("uhv0pth_db", "uhv0pth_db", 200, 0.0, 10.);
1076
1077       //v2pt_db
1078       uhv2pt_db = new TH1D("uhv2pt_db", "uhv2pt_db", 200, 0.0, 10.);
1079       uhv2ptj_db = new TH1D("uhv2ptj_db", "uhv2ptj_db", 200, 0.0, 10.);
1080       uhv2pth_db = new TH1D("uhv2pth_db", "uhv2pth_db", 200, 0.0, 10.);
1081
1082       uhv3pt_db = new TH1D("uhv3pt_db", "uhv3pt_db", 200, 0.0, 10.);
1083       uhv4pt_db = new TH1D("uhv4pt_db", "uhv4pt_db", 200, 0.0, 10.);
1084       uhv5pt_db = new TH1D("uhv5pt_db", "uhv5pt_db", 200, 0.0, 10.);
1085       uhv6pt_db = new TH1D("uhv6pt_db", "uhv6pt_db", 200, 0.0, 10.);
1086
1087       //v0eta
1088       uhv0eta = new TH1D("uhv0eta", "uhv0eta", nintV2eta, v2etaBins);
1089       uhv0etaj = new TH1D("uhv0etaj", "uhv0etaj", nintV2eta, v2etaBins);
1090       uhv0etah = new TH1D("uhv0etah", "uhv0etah", nintV2eta, v2etaBins);
1091
1092       //v2eta
1093       uhv2eta = new TH1D("uhv2eta", "uhv2eta", nintV2eta, v2etaBins);
1094       uhv2etaj = new TH1D("uhv2etaj", "uhv2etaj", nintV2eta, v2etaBins);
1095       uhv2etah = new TH1D("uhv2etah", "uhv2etah", nintV2eta, v2etaBins);
1096
1097       uhv3eta = new TH1D("uhv3eta", "uhv3eta", nintV2eta, v2etaBins);
1098       uhv4eta = new TH1D("uhv4eta", "uhv4eta", nintV2eta, v2etaBins);
1099       uhv5eta = new TH1D("uhv5eta", "uhv5eta", nintV2eta, v2etaBins);
1100       uhv6eta = new TH1D("uhv6eta", "uhv6eta", nintV2eta, v2etaBins);
1101
1102       //v0eta_db
1103       uhv0eta_db = new TH1D("uhv0eta_db", "uhv0eta_db", 200, -5, 5.);
1104       uhv0etaj_db = new TH1D("uhv0etaj_db", "uhv0etaj_db", 200, -5, 5.);
1105       uhv0etah_db = new TH1D("uhv0etah_db", "uhv0etah_db", 200, -5, 5.);
1106
1107       //v2eta_db
1108       uhv2eta_db = new TH1D("uhv2eta_db", "uhv2eta_db", 200, -5, 5.);
1109       uhv2etaj_db = new TH1D("uhv2etaj_db", "uhv2etaj_db", 200, -5, 5.);
1110       uhv2etah_db = new TH1D("uhv2etah_db", "uhv2etah_db", 200, -5, 5.);
1111
1112       uhv3eta_db = new TH1D("uhv3eta_db", "uhv3eta_db", 200, -5, 5.);
1113       uhv4eta_db = new TH1D("uhv4eta_db", "uhv4eta_db", 200, -5, 5.);
1114       uhv5eta_db = new TH1D("uhv5eta_db", "uhv5eta_db", 200, -5, 5.);
1115       uhv6eta_db = new TH1D("uhv6eta_db", "uhv6eta_db", 200, -5, 5.);
1116     }
1117
1118     dhphi = new TH1D("dhphi", "dhphi", 1000, -3.14159265358979, 3.14159265358979);
1119
1120     dhpdg = new TH1D("dhpdg", "dhpdg", 20000001, -10000000.5, 10000000.5);
1121     dhpdg_st = new TH2D("dhpdg_st", "dhpdg_st", 1001, -500.5, 500.5, 3, -0.5, 3.5);
1122
1123     dhet_sum = new TH1D("dhet_sum", "dhet_sum", 300, 0., 100000.);
1124     dhet_barrel_sum = new TH1D("dhet_barrel_sum", "dhet_barrel_sum", 500, 0., 100000.);
1125     dhe_sum = new TH1D("dhe_sum", "dhe_sum", 800, 0., 1000000.);
1126     dhe_barrel_sum = new TH1D("dhe_barrel_sum", "dhe_barrel_sum", 300, 0., 100000.);
1127     dheta = new TH1D("dheta", "dheta", 1000, -10., 10.);
1128     dhpt = new TH1D("dhpt", "dhpt", 1000, 0.0000000000001, 200.);
1129
1130     //charged
1131     dhphi_cha = new TH1D("dhphi_cha", "dhphi_cha", 1000, -3.14159265358979, 3.14159265358979);
1132     dhet_cha_sum = new TH1D("dhet_cha_sum", "dhet_cha_sum", 200, 0., 20000.);
1133     dhet_cha_barrel_sum = new TH1D("dhet_cha_barrel_sum", "dhet_cha_barrel_sum", 300, 0., 10000.);
1134     dheta_cha = new TH1D("dheta_cha", "dheta_cha", 1000, -10., 10.);
1135     dhpt_cha = new TH1D("dhpt_cha", "dhpt_cha", 1000, 0.0000000000001, 100.);
1136
1137     dhv2pt_cha = new TH1D("dhv2pt_cha", "dhv2pt_cha", 200, 0.0, 10.);
1138     dhv0pt_cha = new TH1D("dhv0pt_cha", "dhv0pt_cha", 200, 0.0, 10.);
1139     dhv2eta_cha = new TH1D("dhv2eta_cha", "dhv2eta_cha", 200, -5, 5.);
1140     dhv0eta_cha = new TH1D("dhv0eta_cha", "dhv0eta_cha", 200, -5, 5.);
1141
1142     //charged hadrons
1143     dhphi_ch = new TH1D("dhphi_ch", "dhphi_ch", 1000, -3.14159265358979, 3.14159265358979);
1144     dhet_ch_sum = new TH1D("dhet_ch_sum", "dhet_ch_sum", 200, 0., 20000.);
1145     dhet_ch_barrel_sum = new TH1D("dhet_ch_barrel_sum", "dhet_ch_barrel_sum", 300, 0., 10000.);
1146     dhe_ch_sum = new TH1D("dhe_ch_sum", "dhe_ch_sum", 400, 0., 500000.);
1147     dhe_ch_barrel_sum = new TH1D("dhe_ch_barrel_sum", "dhe_ch_barrel_sum", 150, 0., 10000.);
1148     dheta_ch = new TH1D("dheta_ch", "dheta_ch", 1000, -10., 10.);
1149     dhpt_ch = new TH1D("dhpt_ch", "dhpt_ch", 1000, 0.0000000000001, 100.);
1150
1151     dhv2pt_ch = new TH1D("dhv2pt_ch", "dhv2pt_ch", 200, 0.0, 10.);
1152     dhv0pt_ch = new TH1D("dhv0pt_ch", "dhv0pt_ch", 200, 0.0, 10.);
1153     dhv2eta_ch = new TH1D("dhv2eta_ch", "dhv2eta_ch", 200, -5, 5.);
1154     dhv0eta_ch = new TH1D("dhv0eta_ch", "dhv0eta_ch", 200, -5, 5.);
1155
1156     //ph
1157     dhet_ph_sum = new TH1D("dhet_ph_sum", "dhet_ph_sum", 150, 0., 8000.);
1158     dhet_ph_barrel_sum = new TH1D("dhet_ph_barrel_sum", "dhet_ph_barrel_sum", 100, 0., 5000.);
1159     dhe_ph_sum = new TH1D("dhe_ph_sum", "dhe_ph_sum", 1000, 0., 200000.);
1160     dhe_ph_barrel_sum = new TH1D("dhe_ph_barrel_sum", "dhe_ph_barrel_sum", 100, 0., 5000.);
1161     dheta_ph = new TH1D("dheta_ph", "dheta_ph", 1000, -10., 10.);
1162     dhpt_ph = new TH1D("dhpt_ph", "dhpt_ph", 1000, 0.0000000000001, 100.);
1163
1164     //n
1165     dhet_n_sum = new TH1D("dhet_n_sum", "dhet_n_sum", 150, 0., 3000.);
1166     dhet_n_barrel_sum = new TH1D("dhet_n_barrel_sum", "dhet_n_barrel_sum", 100, 0., 1100.);
1167     dhe_n_sum = new TH1D("dhe_n_sum", "dhe_n_sum", 600, 0., 200000.);
1168     dhe_n_barrel_sum = new TH1D("dhe_n_barrel_sum", "dhe_n_barrel_sum", 100, 0., 1100.);
1169     dheta_n = new TH1D("dheta_n", "dheta_n", 1000, -10., 10.);
1170     dhpt_n = new TH1D("dhpt_n", "dhpt_n", 1000, 0.0000000000001, 100.);
1171
1172     //p
1173     dhet_p_sum = new TH1D("dhet_p_sum", "dhet_p_sum", 150, 0., 3000.);
1174     dhet_p_barrel_sum = new TH1D("dhet_p_barrel_sum", "dhet_p_barrel_sum", 100, 0., 1100.);
1175     dheta_p = new TH1D("dheta_p", "dheta_p", 1000, -10., 10.);
1176     dhpt_p = new TH1D("dhpt_p", "dhpt_p", 1000, 0.0000000000001, 100.);
1177
1178     //pi
1179     dhet_pi_sum = new TH1D("dhet_pi_sum", "dhet_pi_sum", 300, 6000., 18000.);
1180     dhet_pi_barrel_sum = new TH1D("dhet_pi_barrel_sum", "dhet_pi_barrel_sum", 300, 1000., 7000.);
1181     dheta_pi = new TH1D("dheta_pi", "dheta_pi", 1000, -10., 10.);
1182     dhpt_pi = new TH1D("dhpt_pi", "dhpt_pi", 1000, 0.0000000000001, 100.);
1183
1184     //K
1185     dhet_K_sum = new TH1D("dhet_K_sum", "dhet_K_sum", 150, 1500., 4500.);
1186     dhet_K_barrel_sum = new TH1D("dhet_K_barrel_sum", "dhet_K_barrel_sum", 100, 500., 1600.);
1187     dheta_K = new TH1D("dheta_K", "dheta_K", 1000, -10., 10.);
1188     dhpt_K = new TH1D("dhpt_K", "dhpt_K", 1000, 0.0000000000001, 100.);
1189
1190     hNev = new TH1D("hNev", "hNev", 1, 0., 2.);
1191   }
1192
1193   if (doAnalysis_) {
1194     usesResource(TFileService::kSharedResource);
1195     edm::Service<TFileService> f;
1196     nt = f->make<TNtuple>("nt", "Mixing Analysis", "mix:np:src:sig");
1197     hydjetTree_ = f->make<TTree>("hi", "Tree of Hydjet2 Events");
1198     hydjetTree_->Branch("event", &hev_.event, "event/I");
1199     hydjetTree_->Branch("b", &hev_.b, "b/F");
1200     hydjetTree_->Branch("npart", &hev_.npart, "npart/F");
1201     hydjetTree_->Branch("ncoll", &hev_.ncoll, "ncoll/F");
1202     hydjetTree_->Branch("nhard", &hev_.nhard, "nhard/F");
1203     hydjetTree_->Branch("phi0", &hev_.phi0, "phi0/F");
1204     hydjetTree_->Branch("scale", &hev_.scale, "scale/F");
1205     hydjetTree_->Branch("n", hev_.n, "n[3]/I");
1206     hydjetTree_->Branch("ptav", hev_.ptav, "ptav[3]/F");
1207     if (doParticles_) {
1208       hydjetTree_->Branch("mult", &hev_.mult, "mult/I");
1209       hydjetTree_->Branch("px", hev_.px, "px[mult]/F");
1210       hydjetTree_->Branch("py", hev_.py, "py[mult]/F");
1211       hydjetTree_->Branch("pz", hev_.pz, "pz[mult]/F");
1212       hydjetTree_->Branch("e", hev_.e, "e[mult]/F");
1213       hydjetTree_->Branch("pseudoRapidity", hev_.pseudoRapidity, "pseudoRapidity[mult]/F");
1214       hydjetTree_->Branch("pt", hev_.pt, "pt[mult]/F");
1215       hydjetTree_->Branch("eta", hev_.eta, "eta[mult]/F");
1216       hydjetTree_->Branch("phi", hev_.phi, "phi[mult]/F");
1217       hydjetTree_->Branch("pdg", hev_.pdg, "pdg[mult]/I");
1218       hydjetTree_->Branch("chg", hev_.chg, "chg[mult]/I");
1219
1220       hydjetTree_->Branch("vx", &hev_.vx, "vx/F");
1221       hydjetTree_->Branch("vy", &hev_.vy, "vy/F");
1222       hydjetTree_->Branch("vz", &hev_.vz, "vz/F");
1223       hydjetTree_->Branch("vr", &hev_.vr, "vr/F");
1224     }
1225   }
1226 }
1227 // ------------ method called once each job just after ending the event loop  ------------
1228 void Hydjet2Analyzer::endJob() {
1229   if (doHistos_) {
1230     dhphi->Write();
1231     dhpdg->Write();
1232     dhpdg_st->Write();
1233
1234     dhet_sum->Write();
1235     dhet_barrel_sum->Write();
1236     dhe_sum->Write();
1237     dhe_barrel_sum->Write();
1238     dhpt->Write();
1239     dheta->Write();
1240
1241     dhphi_ch->Write();
1242     dhet_ch_sum->Write();
1243     dhet_ch_barrel_sum->Write();
1244     dhe_ch_sum->Write();
1245     dhe_ch_barrel_sum->Write();
1246     dhpt_ch->Write();
1247     dheta_ch->Write();
1248     dhv0pt_ch->Write();
1249     dhv2pt_ch->Write();
1250     dhv0eta_ch->Write();
1251     dhv2eta_ch->Write();
1252
1253     dhphi_cha->Write();
1254     dhet_cha_sum->Write();
1255     dhet_cha_barrel_sum->Write();
1256     dhpt_cha->Write();
1257     dheta_cha->Write();
1258     dhv0pt_cha->Write();
1259     dhv2pt_cha->Write();
1260     dhv0eta_cha->Write();
1261     dhv2eta_cha->Write();
1262
1263     dhet_ph_sum->Write();
1264     dhet_ph_barrel_sum->Write();
1265     dhe_ph_sum->Write();
1266     dhe_ph_barrel_sum->Write();
1267     dhpt_ph->Write();
1268     dheta_ph->Write();
1269
1270     dhet_n_sum->Write();
1271     dhet_n_barrel_sum->Write();
1272     dhe_n_sum->Write();
1273     dhe_n_barrel_sum->Write();
1274     dhpt_n->Write();
1275     dheta_n->Write();
1276
1277     dhet_p_sum->Write();
1278     dhet_p_barrel_sum->Write();
1279     dhpt_p->Write();
1280     dheta_p->Write();
1281
1282     dhet_pi_sum->Write();
1283     dhet_pi_barrel_sum->Write();
1284     dhpt_pi->Write();
1285     dheta_pi->Write();
1286
1287     dhet_K_sum->Write();
1288     dhet_K_barrel_sum->Write();
1289     dhpt_K->Write();
1290     dheta_K->Write();
1291
1292     hNev->Write();
1293     if (userHistos_) {
1294       uhpt->Write();
1295       uhpth->Write();
1296       uhptj->Write();
1297
1298       uhpt_db->Write();
1299       uhpth_db->Write();
1300       uhptj_db->Write();
1301
1302       uhNpart->Write();
1303       uhNparth->Write();
1304       uhNpartj->Write();
1305
1306       uhNpart_db->Write();
1307       uhNparth_db->Write();
1308       uhNpartj_db->Write();
1309
1310       uhPtNpart->Write();
1311       uhPtNparth->Write();
1312       uhPtNpartj->Write();
1313
1314       uhPtNpart_db->Write();
1315       uhPtNparth_db->Write();
1316       uhPtNpartj_db->Write();
1317
1318       uhv2Npart->Write();
1319       uhv2Nparth->Write();
1320       uhv2Npartj->Write();
1321
1322       uhv2Npart_db->Write();
1323       uhv2Nparth_db->Write();
1324       uhv2Npartj_db->Write();
1325
1326       uheta->Write();
1327       uhetah->Write();
1328       uhetaj->Write();
1329       uhphi->Write();
1330       uhphih->Write();
1331       uhphij->Write();
1332
1333       uhv0eta->Write();
1334       uhv0etah->Write();
1335       uhv0etaj->Write();
1336
1337       uhv0eta_db->Write();
1338       uhv0etah_db->Write();
1339       uhv0etaj_db->Write();
1340
1341       uhv0pt->Write();
1342       uhv0pth->Write();
1343       uhv0ptj->Write();
1344
1345       uhv0pt_db->Write();
1346       uhv0pth_db->Write();
1347       uhv0ptj_db->Write();
1348
1349       uhv2eta->Write();
1350       uhv2etah->Write();
1351       uhv2etaj->Write();
1352
1353       uhv2eta_db->Write();
1354       uhv2etah_db->Write();
1355       uhv2etaj_db->Write();
1356
1357       uhv2pt->Write();
1358       uhv2pth->Write();
1359       uhv2ptj->Write();
1360
1361       uhv2pt_db->Write();
1362       uhv2pth_db->Write();
1363       uhv2ptj_db->Write();
1364
1365       uhv3eta->Write();
1366       uhv4eta->Write();
1367       uhv5eta->Write();
1368       uhv6eta->Write();
1369
1370       uhv3eta_db->Write();
1371       uhv4eta_db->Write();
1372       uhv5eta_db->Write();
1373       uhv6eta_db->Write();
1374
1375       uhv3pt->Write();
1376       uhv4pt->Write();
1377       uhv5pt->Write();
1378       uhv6pt->Write();
1379
1380       uhv3pt_db->Write();
1381       uhv4pt_db->Write();
1382       uhv5pt_db->Write();
1383       uhv6pt_db->Write();
1384     }
1385   }
1386 }
1387 //define this as a plug-in
1388 DEFINE_FWK_MODULE(Hydjet2Analyzer);