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