File indexing completed on 2022-03-02 04:16:31
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include <cmath>
0010 #include <fstream>
0011 #include <iostream>
0012 #include <memory>
0013 #include <string>
0014
0015
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
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;
0056 }
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
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;
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
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
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
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
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
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
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
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 }
0406 }
0407 }
0408 Hydjet2Analyzer::~Hydjet2Analyzer() {}
0409
0410
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
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
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
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 {
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();
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
0534
0535
0536
0537
0538
0539
0540
0541
0542
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
0573 if ((*it)->status() < 10) {
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) {
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 }
0594 }
0595
0596
0597 if (userHistos_ && ((uStatus_ == 3) || (((*it)->status() < 10) && (uStatus_ == 1)) ||
0598 (((*it)->status() > 10) && (uStatus_ == 2)))) {
0599
0600
0601 if (std::abs(pdg_id) == uPDG_1 || std::abs(pdg_id) == uPDG_2 || std::abs(pdg_id) == uPDG_3) {
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_) {
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)) {
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)) {
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 }
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)) {
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)) {
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 }
0712
0713 }
0714
0715 }
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 {
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
0798 if (doHistos_) {
0799 hNev->Fill(1., 1);
0800 }
0801 }
0802
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
1093 DEFINE_FWK_MODULE(Hydjet2Analyzer);