JetAnalysis/src/HeavyIonJetAnalyzer.cc

// -*- C++ -*-
//
// Package:    HeavyIonJetAnalyzer
// Class:      HeavyIonJetAnalyzer
// 
/**\class HeavyIonJetAnalyzer HeavyIonJetAnalyzer.cc yetkin/HeavyIonJetAnalyzer/src/HeavyIonJetAnalyzer.cc

 Description: <one line class summary>

 Implementation:
     <Notes on implementation>
*/
//
// Original Author:  Yetkin Yilmaz
//         Created:  Tue Dec 18 09:44:41 EST 2007
// $Id: HeavyIonJetAnalyzer.cc,v 1.7 2009/06/17 16:08:19 yilmaz Exp $
//
//


// system include files
#include <memory>
#include <vector>
#include <iostream>
#include <string>
#include <fstream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/EventSetup.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/InputTag.h"
#include "FWCore/Framework/interface/ESHandle.h"

#include "FWCore/ServiceRegistry/interface/Service.h"

#include "PhysicsTools/UtilAlgos/interface/TFileService.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"

#include "DataFormats/JetReco/interface/JetCollection.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/JetReco/interface/GenJetCollection.h"
#include "DataFormats/JetReco/interface/GenJet.h"


#include "HepMC/GenEvent.h"
#include "HepMC/HeavyIon.h"

#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"

// root include file
#include "TFile.h"
#include "TNtuple.h"

#include "CmsHi/JetAnalysis/macros/CmsHiFunctions.h"

using namespace std;

//static const int PI = 3.14159265358979;

static const int MAXPARTICLES = 5000000;
static const int MAXJETS = 500000;
static const int MAXCONS = 200;

static const int MAXHITS = 50000;
static const int MAXVTX = 1000;
static const int ETABINS = 3; // Fix also in branch string

//
// class decleration
//

struct JetAxis{
  const reco::Jet* axis_;
  bool operator() (const reco::Candidate * cand1,const reco::Candidate * cand2){ 
    double dr1 = deltaR(axis_->eta(),axis_->phi(),cand1->eta(),cand1->phi());
    double dr2 = deltaR(axis_->eta(),axis_->phi(),cand2->eta(),cand2->phi());

    return dr1 < dr2;
  }
};

struct MyJet{

  float etjet;
  float etajet;
  float phijet;
  float area;
  int ncons;

  float r20;
  float r50;
  float r90;

  float e01;
  float e02;
  float e03;
  float e04;
  float e05;

  float et[MAXCONS];
  float eta[MAXCONS];
  float phi[MAXCONS];

  float dr[MAXCONS];

  int event;
  float b;

};


struct HydjetEvent{

   int event;
   float b;
   float npart;
   float ncoll;
   float nhard;
   float phi0;

   int n[ETABINS];
   float ptav[ETABINS];

   int np;
   float par_pt[MAXPARTICLES];
   float par_eta[MAXPARTICLES];
   float par_phi[MAXPARTICLES];
   int pdg[MAXPARTICLES];
   int chg[MAXPARTICLES];

   int algos;
   int njet;

   float et[MAXJETS];
   float eta[MAXJETS];
   float phi[MAXJETS];
   float area[MAXJETS];

  int ncons[MAXJETS];

  float r20[MAXJETS];
  float r50[MAXJETS];
  float r90[MAXJETS];

  float e01[MAXJETS];
  float e02[MAXJETS];
  float e03[MAXJETS];
  float e04[MAXJETS];
  float e05[MAXJETS];
 
   float vx;
   float vy;
   float vz;
   float vr;

};

class HeavyIonJetAnalyzer : public edm::EDAnalyzer {
   public:
      explicit HeavyIonJetAnalyzer(const edm::ParameterSet&);
      ~HeavyIonJetAnalyzer();
    bool sortDeltaR(const reco::Candidate * cand1,const reco::Candidate * cand2);

   private:
      virtual void beginJob(const edm::EventSetup&) ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

      // ----------member data ---------------------------

   std::ofstream out_b;
   std::string fBFileName;

   std::ofstream out_n;
   std::string fNFileName;

   std::ofstream out_m;
   std::string fMFileName;

  
   TTree* hydjetTree_;
  TTree* jetTree_;

   HydjetEvent hev_;
  MyJet yet_;

   TNtuple *ntpart;

   std::string output;           // Output filename
 
   bool doAnalysis_;
   bool doParticles_;
   bool printLists_;
   bool doCF_;
   bool doVertices_;
   double etaMax_;
   double ptMin_;

   edm::InputTag jetSrc_;
   string hepmcSrc_;
   edm::ESHandle < ParticleDataTable > pdt;
   edm::Service<TFileService> f;

};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
HeavyIonJetAnalyzer::HeavyIonJetAnalyzer(const edm::ParameterSet& iConfig)
{
   //now do what ever initialization is needed
   fBFileName = iConfig.getUntrackedParameter<std::string>("output_b", "b_values.txt");
   fNFileName = iConfig.getUntrackedParameter<std::string>("output_n", "n_values.txt");
   fMFileName = iConfig.getUntrackedParameter<std::string>("output_m", "m_values.txt");
   doAnalysis_ = iConfig.getUntrackedParameter<bool>("doAnalysis", true);
   doParticles_ = iConfig.getUntrackedParameter<bool>("doParticles", true);
   printLists_ = iConfig.getUntrackedParameter<bool>("printLists", false);
   doCF_ = iConfig.getUntrackedParameter<bool>("doMixed", false);
   doVertices_ = iConfig.getUntrackedParameter<bool>("doVertices", false);
   jetSrc_ = iConfig.getParameter<edm::InputTag>("jetSrc");
   hepmcSrc_ = iConfig.getUntrackedParameter<string>("hepmcSrc","generator");

   etaMax_ = iConfig.getUntrackedParameter<double>("etaMax", 2);
   ptMin_ = iConfig.getUntrackedParameter<double>("ptMin", 0);

   // Output

}


HeavyIonJetAnalyzer::~HeavyIonJetAnalyzer()
{
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}

//
// member functions
//

// ------------ method called to for each event  ------------
void
HeavyIonJetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{

  cout<<"Analyze"<<endl;

   using namespace edm;
   using namespace HepMC;
  
   hev_.event = iEvent.id().event();
   for(int ieta = 0; ieta < ETABINS; ++ieta){
      hev_.n[ieta] = 0;
      hev_.ptav[ieta] = 0;
   }

   hev_.njet = 0;
   hev_.np = 0;

   double b = -1;
   int npart = -1;
   int ncoll = -1;
   int nhard = -1;
   double phi0 = -99;

   double vx = -99;
   double vy = -99;
   double vz = -99;
   double vr = -99;
   const GenEvent* evt;
   
   if(doCF_){

     Handle<CrossingFrame<HepMCProduct> > cf;
     iEvent.getByLabel(InputTag("mix","source"),cf);

     /*

     MixCollection<HepMCProduct> mix(cf.product());

     int mixsize = mix.size();

     cout<<"Mix Collection Size: "<<mixsize<<endl;
     evt = mix.getObject(mixsize-1).GetEvent();

     MixCollection<HepMCProduct>::iterator begin = mix.begin();
     MixCollection<HepMCProduct>::iterator end = mix.end();
     
     for(MixCollection<HepMCProduct>::iterator mixit = begin; mixit != end; ++mixit){

       const GenEvent* subevt = (*mixit).GetEvent();
       int all = subevt->particles_size();
       HepMC::GenEvent::particle_const_iterator begin = subevt->particles_begin();
       HepMC::GenEvent::particle_const_iterator end = subevt->particles_end();
       for(HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it){
         if((*it)->status() == 1){
           float pdg_id = (*it)->pdg_id();
           float eta = (*it)->momentum().eta();
           float pt = (*it)->momentum().perp();
           const ParticleData * part = pdt->particle(pdg_id );
           float charge = part->charge();
           }
         }
       }
       
     }     

     */

   }else{
      
      Handle<HepMCProduct> mc;
      iEvent.getByLabel(hepmcSrc_,mc);
      evt = mc->GetEvent();

      if(doParticles_){
      int all = evt->particles_size();
      HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
      HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
      for(HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it){
        if((*it)->status() == 1){
           int pdg_id = (*it)->pdg_id();
           float eta = (*it)->momentum().eta();
           float phi = (*it)->momentum().phi();
           float pt = (*it)->momentum().perp();
          const ParticleData * part = pdt->particle(pdg_id );
          int charge = part->charge();

          hev_.par_pt[hev_.np] = pt;
          hev_.par_eta[hev_.np] = eta;
          hev_.par_phi[hev_.np] = phi;
          hev_.pdg[hev_.np] = pdg_id;
          hev_.chg[hev_.np] = charge;

          eta = fabs(eta);
          int etabin = 0;
          if(eta > 0.5) etabin = 1; 
          if(eta > 1.) etabin = 2;
          if(eta < 2.){
             hev_.ptav[etabin] += pt;
             ++(hev_.n[etabin]);
          }
          ++(hev_.np);
        }
      }
      }

   }

   const HeavyIon* hi = evt->heavy_ion();
   if(hi){
      b = hi->impact_parameter();
      npart = hi->Npart_proj()+hi->Npart_targ();
      ncoll = hi->Ncoll();
      nhard = hi->Ncoll_hard();
      phi0 = hi->event_plane_angle();

      if(printLists_){
         out_b<<b<<endl;
         out_n<<npart<<endl;
      }

   }
   
   //   edm::Handle<reco::GenJetCollection> genjets;
   edm::Handle<reco::JetView> genjets;
   //   edm::Handle<vector<reco::Jet> > genjets;
   iEvent.getByLabel(jetSrc_,genjets);
   for(int ijet = 0; ijet < genjets->size(); ++ijet){

      const reco::Jet* jet = &((*genjets)[ijet]);
      //         const reco::GenJet* jet = dynamic_cast<const reco::GenJet*>(&((*genjets)[ijet]));

      cout<<"Jet Quantities : "<<jet->et()<<" "<<jet->eta()<<" "<<jet->phi()<<" "<<jet->jetArea()<<endl;

      double phi = jet->phi()-phi0;
      if(phi > PI) phi = phi - 2*PI;
      if(phi < -PI) phi = phi + 2*PI;
      
      hev_.et[hev_.njet] = jet->et();
      hev_.eta[hev_.njet] = jet->eta();
      hev_.phi[hev_.njet] = phi;

      hev_.area[hev_.njet] = jet->jetArea();

      // Constituents
      std::vector<const reco::Candidate*> members = jet->getJetConstituentsQuick();

      int nm = members.size();

      hev_.ncons[hev_.njet] = nm;
      yet_.ncons = nm;

      cout<<"Jet has "<<nm<<" constituents."<<endl;
      
      JetAxis jaxis;
      jaxis.axis_ = jet;
      
      //      cout<<"Jet axis is at : "<<jaxis.axis_->eta()<<"  "<<jaxis.axis_->phi()<<endl;

      sort(members.begin(),members.end(),jaxis);

      double sum = 0;
      double et = jet->et();
      double r20 = -99;
      double r50 = -99;
      double r90 = -99;

      double e0[6] = {0,0,0,0,0,0};

      for(int im = 0; im < nm; ++im){
        const reco::Candidate* candi = members[im];
        double dr = deltaR(jet->eta(),jet->phi(),candi->eta(),candi->phi());
        cout<<"Constiutuent's Delta R = "<<dr<<endl;

        double phicon = candi->phi()-phi0;
        if(phicon > PI) phicon = phicon - 2*PI;
        if(phicon < -PI) phicon = phicon + 2*PI;


        yet_.et[im] = candi->et();
        yet_.eta[im] = candi->eta();
        yet_.phi[im] = phicon;  
        yet_.dr[im] = dr;

        for(int ir = 1; ir <6; ++ir){
          if(dr < 0.1*ir){
            e0[ir] += candi->et();        
          }
        }
        
        sum += candi->et();
        
        if(sum > et*0.2 && r20 == -99) r20 = dr;
        if(sum > et*0.5 && r50 == -99) r50 = dr;
        if(sum > et*0.9 && r90 == -99) r90 = dr;
        
      }

      hev_.r20[hev_.njet] = r20;
      hev_.r50[hev_.njet] = r50;
      hev_.r90[hev_.njet] = r90;

      hev_.e01[hev_.njet] = e0[1]/et;
      hev_.e02[hev_.njet] = e0[2]/et;
      hev_.e03[hev_.njet] = e0[3]/et;
      hev_.e04[hev_.njet] = e0[4]/et;
      hev_.e05[hev_.njet] = e0[5]/et;

      yet_.e01 = e0[1]/et;
      yet_.e02 = e0[2]/et;
      yet_.e03 = e0[3]/et;
      yet_.e04 = e0[4]/et;
      yet_.e05 = e0[5]/et;

      yet_.etjet = jet->et();
      yet_.etajet = jet->eta();
      yet_.phijet = phi;
      yet_.area = jet->jetArea();

      yet_.event = hev_.event;
      yet_.b = hev_.b;

      jetTree_->Fill();

      ++(hev_.njet);
      
   }
   
   if(doVertices_){   

     HepMC::GenVertex* genvtx = evt->signal_process_vertex();
     if(!genvtx){
       cout<<"No Signal Process Vertex!"<<endl;
       HepMC::GenEvent::particle_const_iterator pit=evt->particles_begin();
       HepMC::GenEvent::particle_const_iterator ptend=evt->particles_end();
       while(!genvtx || ( genvtx->particles_in_size() == 1 && pit != ptend ) ){
         if(!genvtx) cout<<"No Gen Vertex!"<<endl;
         ++pit;
         if(pit == ptend) cout<<"End reached!"<<endl;
         genvtx = (*pit)->production_vertex();
       }
     }

   vx = genvtx->position().x()/10;
   vy = genvtx->position().y()/10;
   vz = genvtx->position().z()/10;
   vr = genvtx->position().rho()/10;

   }

   hev_.b = b;
   hev_.npart = npart;
   hev_.ncoll = ncoll;
   hev_.nhard = nhard;
   hev_.phi0 = phi0;

   hev_.vx = vx;
   hev_.vy = vy;
   hev_.vz = vz;
   hev_.vr = vr;

   hydjetTree_->Fill();

}


// ------------ method called once each job just before starting event loop  ------------
void 
HeavyIonJetAnalyzer::beginJob(const edm::EventSetup& iSetup)
{
  cout<<"Begin"<<endl;

   iSetup.getData(pdt);

   if(printLists_){
      out_b.open(fBFileName.c_str());
      if(out_b.good() == false)
         throw cms::Exception("BadFile") << "Can\'t open file " << fBFileName;
      out_n.open(fNFileName.c_str());
      if(out_n.good() == false)
         throw cms::Exception("BadFile") << "Can\'t open file " << fNFileName;
      out_m.open(fMFileName.c_str());
      if(out_m.good() == false)
         throw cms::Exception("BadFile") << "Can\'t open file " << fMFileName;
   }   
   
   if(doAnalysis_){
      hydjetTree_ = f->make<TTree>("hi","Tree of Hydjet Events");
      jetTree_ = f->make<TTree>("jet","Tree of Jet Constituents");

      hydjetTree_->Branch("event",&hev_.event,"event/I");
      hydjetTree_->Branch("b",&hev_.b,"b/F");
      hydjetTree_->Branch("npart",&hev_.npart,"npart/F");
      hydjetTree_->Branch("ncoll",&hev_.ncoll,"ncoll/F");
      hydjetTree_->Branch("nhard",&hev_.nhard,"nhard/F");
      hydjetTree_->Branch("phi0",&hev_.phi0,"phi0/F");

      if(doParticles_){
         hydjetTree_->Branch("n",hev_.n,"n[3]/I");
         hydjetTree_->Branch("ptav",hev_.ptav,"ptav[3]/F");
         hydjetTree_->Branch("np",&hev_.np,"np/I");
         hydjetTree_->Branch("par_pt",hev_.par_pt,"par_pt[np]/F");
         hydjetTree_->Branch("par_eta",hev_.par_eta,"par_eta[np]/F");
         hydjetTree_->Branch("par_phi",hev_.par_phi,"par_phi[np]/F");
         hydjetTree_->Branch("pdg",hev_.pdg,"pdg[np]/I");
         hydjetTree_->Branch("chg",hev_.chg,"chg[np]/I");
      }

      hydjetTree_->Branch("vx",&hev_.vx,"vx/F");
      hydjetTree_->Branch("vy",&hev_.vy,"vy/F");
      hydjetTree_->Branch("vz",&hev_.vz,"vz/F");
      hydjetTree_->Branch("vr",&hev_.vr,"vr/F");

      hydjetTree_->Branch("njet",&hev_.njet,"njet/I");
      hydjetTree_->Branch("et",hev_.et,"et[njet]/F");
      hydjetTree_->Branch("eta",hev_.eta,"eta[njet]/F");
      hydjetTree_->Branch("phi",hev_.phi,"phi[njet]/F");
      hydjetTree_->Branch("area",hev_.area,"area[njet]/F");

      hydjetTree_->Branch("r20",hev_.r20,"r20[njet]/F");
      hydjetTree_->Branch("r50",hev_.r50,"r50[njet]/F");
      hydjetTree_->Branch("r90",hev_.r90,"r90[njet]/F");

      hydjetTree_->Branch("e01",hev_.e01,"e01[njet]/F");
      hydjetTree_->Branch("e02",hev_.e02,"e02[njet]/F");
      hydjetTree_->Branch("e03",hev_.e03,"e03[njet]/F");
      hydjetTree_->Branch("e04",hev_.e04,"e04[njet]/F");
      hydjetTree_->Branch("e05",hev_.e05,"e05[njet]/F");

      hydjetTree_->Branch("ncons",hev_.ncons,"ncons[njet]/I");

      jetTree_->Branch("event",&hev_.event,"event/I");
      jetTree_->Branch("b",&hev_.b,"b/F");

      jetTree_->Branch("ncons",&yet_.ncons,"ncons/I");
      jetTree_->Branch("et",yet_.et,"et[ncons]/F");
      jetTree_->Branch("eta",yet_.eta,"eta[ncons]/F");
      jetTree_->Branch("phi",yet_.phi,"phi[ncons]/F");
      jetTree_->Branch("dr",yet_.dr,"dr[ncons]/F");

      jetTree_->Branch("etjet",&yet_.etjet,"etjet/F");
      jetTree_->Branch("etajet",&yet_.etajet,"etajet/F");
      jetTree_->Branch("phijet",&yet_.phijet,"phijet/F");

      jetTree_->Branch("area",&yet_.area,"area/F");
      jetTree_->Branch("r20",&yet_.r20,"r20/F");
      jetTree_->Branch("r50",&yet_.r50,"r50/F");
      jetTree_->Branch("r90",&yet_.r90,"r90/F");

      jetTree_->Branch("e01",&yet_.e01,"e01/F");
      jetTree_->Branch("e02",&yet_.e02,"e02/F");
      jetTree_->Branch("e03",&yet_.e03,"e03/F");
      jetTree_->Branch("e04",&yet_.e04,"e04/F");
      jetTree_->Branch("e05",&yet_.e05,"e05/F");


   }
  
}

// ------------ method called once each job just after ending the event loop  ------------
void 
HeavyIonJetAnalyzer::endJob() {
}

//define this as a plug-in
DEFINE_FWK_MODULE(HeavyIonJetAnalyzer);
Revision:	1.8
Committed:	Tue Jul 21 14:43:46 2009 UTC (15 years, 9 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.7:	+4 -4 lines
Log Message:	GenJet cleaner for matching
#	User	Rev	Content
1	yilmaz	1.1	// -- C++ --
2			//
3			// Package: HeavyIonJetAnalyzer
4			// Class: HeavyIonJetAnalyzer
5			//
6			/**\class HeavyIonJetAnalyzer HeavyIonJetAnalyzer.cc yetkin/HeavyIonJetAnalyzer/src/HeavyIonJetAnalyzer.cc
7
8			Description: <one line class summary>
9
10			Implementation:
11			<Notes on implementation>
12			*/
13			//
14			// Original Author: Yetkin Yilmaz
15			// Created: Tue Dec 18 09:44:41 EST 2007
16	yilmaz	1.8	// $Id: HeavyIonJetAnalyzer.cc,v 1.7 2009/06/17 16:08:19 yilmaz Exp $
17	yilmaz	1.1	//
18			//
19
20
21			// system include files
22			#include <memory>
23			#include <vector>
24			#include <iostream>
25			#include <string>
26			#include <fstream>
27
28			// user include files
29			#include "FWCore/Framework/interface/Frameworkfwd.h"
30			#include "FWCore/Framework/interface/EDAnalyzer.h"
31
32			#include "FWCore/Framework/interface/Event.h"
33			#include "FWCore/Framework/interface/MakerMacros.h"
34			#include "FWCore/Framework/interface/EventSetup.h"
35
36			#include "FWCore/ParameterSet/interface/ParameterSet.h"
37			#include "FWCore/ParameterSet/interface/InputTag.h"
38			#include "FWCore/Framework/interface/ESHandle.h"
39
40			#include "FWCore/ServiceRegistry/interface/Service.h"
41
42			#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
43
44	yilmaz	1.5	#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
45	yilmaz	1.1	#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
46
47			#include "DataFormats/JetReco/interface/JetCollection.h"
48			#include "DataFormats/JetReco/interface/Jet.h"
49			#include "DataFormats/JetReco/interface/GenJetCollection.h"
50			#include "DataFormats/JetReco/interface/GenJet.h"
51
52
53			#include "HepMC/GenEvent.h"
54			#include "HepMC/HeavyIon.h"
55
56			#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
57
58			// root include file
59			#include "TFile.h"
60			#include "TNtuple.h"
61
62	yilmaz	1.4	#include "CmsHi/JetAnalysis/macros/CmsHiFunctions.h"
63
64	yilmaz	1.1	using namespace std;
65
66	yilmaz	1.4	//static const int PI = 3.14159265358979;
67	yilmaz	1.1
68	yilmaz	1.4	static const int MAXPARTICLES = 5000000;
69			static const int MAXJETS = 500000;
70			static const int MAXCONS = 200;
71
72			static const int MAXHITS = 50000;
73			static const int MAXVTX = 1000;
74			static const int ETABINS = 3; // Fix also in branch string
75	yilmaz	1.1
76			//
77			// class decleration
78			//
79
80	yilmaz	1.4	struct JetAxis{
81			const reco::Jet* axis_;
82			bool operator() (const reco::Candidate * cand1,const reco::Candidate * cand2){
83			double dr1 = deltaR(axis_->eta(),axis_->phi(),cand1->eta(),cand1->phi());
84			double dr2 = deltaR(axis_->eta(),axis_->phi(),cand2->eta(),cand2->phi());
85
86			return dr1 < dr2;
87			}
88			};
89
90			struct MyJet{
91
92			float etjet;
93			float etajet;
94			float phijet;
95			float area;
96			int ncons;
97
98			float r20;
99			float r50;
100			float r90;
101
102			float e01;
103			float e02;
104			float e03;
105			float e04;
106			float e05;
107
108			float et[MAXCONS];
109			float eta[MAXCONS];
110			float phi[MAXCONS];
111
112			float dr[MAXCONS];
113
114			int event;
115			float b;
116
117			};
118
119
120	yilmaz	1.1	struct HydjetEvent{
121
122			int event;
123			float b;
124			float npart;
125			float ncoll;
126			float nhard;
127	yilmaz	1.4	float phi0;
128	yilmaz	1.1
129			int n[ETABINS];
130			float ptav[ETABINS];
131
132			int np;
133			float par_pt[MAXPARTICLES];
134			float par_eta[MAXPARTICLES];
135			float par_phi[MAXPARTICLES];
136			int pdg[MAXPARTICLES];
137			int chg[MAXPARTICLES];
138
139			int algos;
140			int njet;
141
142			float et[MAXJETS];
143			float eta[MAXJETS];
144			float phi[MAXJETS];
145			float area[MAXJETS];
146	yilmaz	1.4
147			int ncons[MAXJETS];
148
149			float r20[MAXJETS];
150			float r50[MAXJETS];
151			float r90[MAXJETS];
152
153			float e01[MAXJETS];
154			float e02[MAXJETS];
155			float e03[MAXJETS];
156			float e04[MAXJETS];
157			float e05[MAXJETS];
158
159	yilmaz	1.1	float vx;
160			float vy;
161			float vz;
162			float vr;
163
164			};
165
166			class HeavyIonJetAnalyzer : public edm::EDAnalyzer {
167			public:
168			explicit HeavyIonJetAnalyzer(const edm::ParameterSet&);
169			~HeavyIonJetAnalyzer();
170	yilmaz	1.4	bool sortDeltaR(const reco::Candidate * cand1,const reco::Candidate * cand2);
171	yilmaz	1.1
172			private:
173			virtual void beginJob(const edm::EventSetup&) ;
174			virtual void analyze(const edm::Event&, const edm::EventSetup&);
175			virtual void endJob() ;
176
177			// ----------member data ---------------------------
178
179			std::ofstream out_b;
180			std::string fBFileName;
181
182			std::ofstream out_n;
183			std::string fNFileName;
184
185			std::ofstream out_m;
186			std::string fMFileName;
187
188
189			TTree* hydjetTree_;
190	yilmaz	1.4	TTree* jetTree_;
191
192	yilmaz	1.1	HydjetEvent hev_;
193	yilmaz	1.4	MyJet yet_;
194	yilmaz	1.1
195			TNtuple *ntpart;
196
197			std::string output; // Output filename
198
199			bool doAnalysis_;
200			bool doParticles_;
201			bool printLists_;
202			bool doCF_;
203	yilmaz	1.4	bool doVertices_;
204	yilmaz	1.1	double etaMax_;
205			double ptMin_;
206
207	yilmaz	1.8	edm::InputTag jetSrc_;
208	yilmaz	1.7	string hepmcSrc_;
209	yilmaz	1.1	edm::ESHandle < ParticleDataTable > pdt;
210			edm::Service<TFileService> f;
211
212			};
213
214			//
215			// constants, enums and typedefs
216			//
217
218			//
219			// static data member definitions
220			//
221
222			//
223			// constructors and destructor
224			//
225			HeavyIonJetAnalyzer::HeavyIonJetAnalyzer(const edm::ParameterSet& iConfig)
226			{
227			//now do what ever initialization is needed
228			fBFileName = iConfig.getUntrackedParameter<std::string>("output_b", "b_values.txt");
229			fNFileName = iConfig.getUntrackedParameter<std::string>("output_n", "n_values.txt");
230			fMFileName = iConfig.getUntrackedParameter<std::string>("output_m", "m_values.txt");
231			doAnalysis_ = iConfig.getUntrackedParameter<bool>("doAnalysis", true);
232			doParticles_ = iConfig.getUntrackedParameter<bool>("doParticles", true);
233			printLists_ = iConfig.getUntrackedParameter<bool>("printLists", false);
234			doCF_ = iConfig.getUntrackedParameter<bool>("doMixed", false);
235	yilmaz	1.4	doVertices_ = iConfig.getUntrackedParameter<bool>("doVertices", false);
236	yilmaz	1.8	jetSrc_ = iConfig.getParameter<edm::InputTag>("jetSrc");
237	yilmaz	1.7	hepmcSrc_ = iConfig.getUntrackedParameter<string>("hepmcSrc","generator");
238	yilmaz	1.1
239			etaMax_ = iConfig.getUntrackedParameter<double>("etaMax", 2);
240			ptMin_ = iConfig.getUntrackedParameter<double>("ptMin", 0);
241
242			// Output
243
244			}
245
246
247			HeavyIonJetAnalyzer::~HeavyIonJetAnalyzer()
248			{
249			// do anything here that needs to be done at desctruction time
250			// (e.g. close files, deallocate resources etc.)
251
252			}
253
254			//
255			// member functions
256			//
257
258			// ------------ method called to for each event ------------
259			void
260			HeavyIonJetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
261			{
262	yilmaz	1.4
263			cout<<"Analyze"<<endl;
264
265	yilmaz	1.1	using namespace edm;
266			using namespace HepMC;
267
268			hev_.event = iEvent.id().event();
269			for(int ieta = 0; ieta < ETABINS; ++ieta){
270			hev_.n[ieta] = 0;
271			hev_.ptav[ieta] = 0;
272			}
273
274			hev_.njet = 0;
275			hev_.np = 0;
276
277			double b = -1;
278			int npart = -1;
279			int ncoll = -1;
280			int nhard = -1;
281	yilmaz	1.4	double phi0 = -99;
282
283	yilmaz	1.1	double vx = -99;
284			double vy = -99;
285			double vz = -99;
286			double vr = -99;
287			const GenEvent* evt;
288	yilmaz	1.2
289	yilmaz	1.1	if(doCF_){
290
291			Handle<CrossingFrame<HepMCProduct> > cf;
292			iEvent.getByLabel(InputTag("mix","source"),cf);
293
294			/*
295
296			MixCollection<HepMCProduct> mix(cf.product());
297
298			int mixsize = mix.size();
299
300			cout<<"Mix Collection Size: "<<mixsize<<endl;
301			evt = mix.getObject(mixsize-1).GetEvent();
302
303			MixCollection<HepMCProduct>::iterator begin = mix.begin();
304			MixCollection<HepMCProduct>::iterator end = mix.end();
305
306			for(MixCollection<HepMCProduct>::iterator mixit = begin; mixit != end; ++mixit){
307
308			const GenEvent* subevt = (*mixit).GetEvent();
309			int all = subevt->particles_size();
310			HepMC::GenEvent::particle_const_iterator begin = subevt->particles_begin();
311			HepMC::GenEvent::particle_const_iterator end = subevt->particles_end();
312			for(HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it){
313			if((*it)->status() == 1){
314			float pdg_id = (*it)->pdg_id();
315			float eta = (*it)->momentum().eta();
316			float pt = (*it)->momentum().perp();
317			const ParticleData * part = pdt->particle(pdg_id );
318			float charge = part->charge();
319			}
320			}
321			}
322
323			}
324
325			*/
326
327			}else{
328
329			Handle<HepMCProduct> mc;
330	yilmaz	1.7	iEvent.getByLabel(hepmcSrc_,mc);
331	yilmaz	1.1	evt = mc->GetEvent();
332
333			if(doParticles_){
334	yilmaz	1.2	int all = evt->particles_size();
335			HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
336			HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
337			for(HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it){
338			if((*it)->status() == 1){
339			int pdg_id = (*it)->pdg_id();
340			float eta = (*it)->momentum().eta();
341			float phi = (*it)->momentum().phi();
342			float pt = (*it)->momentum().perp();
343			const ParticleData * part = pdt->particle(pdg_id );
344			int charge = part->charge();
345
346			hev_.par_pt[hev_.np] = pt;
347			hev_.par_eta[hev_.np] = eta;
348			hev_.par_phi[hev_.np] = phi;
349			hev_.pdg[hev_.np] = pdg_id;
350			hev_.chg[hev_.np] = charge;
351
352			eta = fabs(eta);
353			int etabin = 0;
354			if(eta > 0.5) etabin = 1;
355			if(eta > 1.) etabin = 2;
356			if(eta < 2.){
357			hev_.ptav[etabin] += pt;
358			++(hev_.n[etabin]);
359			}
360			++(hev_.np);
361			}
362			}
363			}
364	yilmaz	1.1
365			}
366
367			const HeavyIon* hi = evt->heavy_ion();
368			if(hi){
369			b = hi->impact_parameter();
370			npart = hi->Npart_proj()+hi->Npart_targ();
371	yilmaz	1.3	ncoll = hi->Ncoll();
372			nhard = hi->Ncoll_hard();
373	yilmaz	1.4	phi0 = hi->event_plane_angle();
374	yilmaz	1.1
375			if(printLists_){
376			out_b<<b<<endl;
377			out_n<<npart<<endl;
378			}
379
380			}
381
382			// edm::Handle<reco::GenJetCollection> genjets;
383			edm::Handle<reco::JetView> genjets;
384			// edm::Handle<vector<reco::Jet> > genjets;
385	yilmaz	1.8	iEvent.getByLabel(jetSrc_,genjets);
386	yilmaz	1.1	for(int ijet = 0; ijet < genjets->size(); ++ijet){
387
388			const reco::Jet* jet = &((*genjets)[ijet]);
389			// const reco::GenJet* jet = dynamic_cast<const reco::GenJet>(&((genjets)[ijet]));
390
391	yilmaz	1.4	cout<<"Jet Quantities : "<<jet->et()<<" "<<jet->eta()<<" "<<jet->phi()<<" "<<jet->jetArea()<<endl;
392
393			double phi = jet->phi()-phi0;
394			if(phi > PI) phi = phi - 2*PI;
395			if(phi < -PI) phi = phi + 2*PI;
396	yilmaz	1.1
397	yilmaz	1.4	hev_.et[hev_.njet] = jet->et();
398	yilmaz	1.1	hev_.eta[hev_.njet] = jet->eta();
399	yilmaz	1.4	hev_.phi[hev_.njet] = phi;
400
401			hev_.area[hev_.njet] = jet->jetArea();
402
403			// Constituents
404			std::vector<const reco::Candidate*> members = jet->getJetConstituentsQuick();
405
406			int nm = members.size();
407
408			hev_.ncons[hev_.njet] = nm;
409			yet_.ncons = nm;
410
411			cout<<"Jet has "<<nm<<" constituents."<<endl;
412
413			JetAxis jaxis;
414			jaxis.axis_ = jet;
415
416			// cout<<"Jet axis is at : "<<jaxis.axis_->eta()<<" "<<jaxis.axis_->phi()<<endl;
417
418			sort(members.begin(),members.end(),jaxis);
419
420			double sum = 0;
421			double et = jet->et();
422			double r20 = -99;
423			double r50 = -99;
424			double r90 = -99;
425
426			double e0[6] = {0,0,0,0,0,0};
427
428			for(int im = 0; im < nm; ++im){
429			const reco::Candidate* candi = members[im];
430			double dr = deltaR(jet->eta(),jet->phi(),candi->eta(),candi->phi());
431			cout<<"Constiutuent's Delta R = "<<dr<<endl;
432
433			double phicon = candi->phi()-phi0;
434			if(phicon > PI) phicon = phicon - 2*PI;
435			if(phicon < -PI) phicon = phicon + 2*PI;
436
437
438			yet_.et[im] = candi->et();
439			yet_.eta[im] = candi->eta();
440			yet_.phi[im] = phicon;
441			yet_.dr[im] = dr;
442
443			for(int ir = 1; ir <6; ++ir){
444			if(dr < 0.1*ir){
445			e0[ir] += candi->et();
446			}
447			}
448
449			sum += candi->et();
450
451			if(sum > et*0.2 && r20 == -99) r20 = dr;
452			if(sum > et*0.5 && r50 == -99) r50 = dr;
453			if(sum > et*0.9 && r90 == -99) r90 = dr;
454
455			}
456
457			hev_.r20[hev_.njet] = r20;
458			hev_.r50[hev_.njet] = r50;
459			hev_.r90[hev_.njet] = r90;
460
461			hev_.e01[hev_.njet] = e0[1]/et;
462			hev_.e02[hev_.njet] = e0[2]/et;
463			hev_.e03[hev_.njet] = e0[3]/et;
464			hev_.e04[hev_.njet] = e0[4]/et;
465			hev_.e05[hev_.njet] = e0[5]/et;
466
467			yet_.e01 = e0[1]/et;
468			yet_.e02 = e0[2]/et;
469			yet_.e03 = e0[3]/et;
470			yet_.e04 = e0[4]/et;
471			yet_.e05 = e0[5]/et;
472
473			yet_.etjet = jet->et();
474			yet_.etajet = jet->eta();
475			yet_.phijet = phi;
476			yet_.area = jet->jetArea();
477
478			yet_.event = hev_.event;
479			yet_.b = hev_.b;
480
481			jetTree_->Fill();
482	yilmaz	1.1
483			++(hev_.njet);
484	yilmaz	1.2
485	yilmaz	1.1	}
486
487	yilmaz	1.4	if(doVertices_){
488	yilmaz	1.6
489			HepMC::GenVertex* genvtx = evt->signal_process_vertex();
490			if(!genvtx){
491			cout<<"No Signal Process Vertex!"<<endl;
492			HepMC::GenEvent::particle_const_iterator pit=evt->particles_begin();
493			HepMC::GenEvent::particle_const_iterator ptend=evt->particles_end();
494			while(!genvtx \|\| ( genvtx->particles_in_size() == 1 && pit != ptend ) ){
495			if(!genvtx) cout<<"No Gen Vertex!"<<endl;
496			++pit;
497			if(pit == ptend) cout<<"End reached!"<<endl;
498			genvtx = (*pit)->production_vertex();
499			}
500			}
501
502			vx = genvtx->position().x()/10;
503			vy = genvtx->position().y()/10;
504			vz = genvtx->position().z()/10;
505			vr = genvtx->position().rho()/10;
506
507	yilmaz	1.4	}
508	yilmaz	1.1
509			hev_.b = b;
510			hev_.npart = npart;
511			hev_.ncoll = ncoll;
512			hev_.nhard = nhard;
513	yilmaz	1.4	hev_.phi0 = phi0;
514
515	yilmaz	1.1	hev_.vx = vx;
516			hev_.vy = vy;
517			hev_.vz = vz;
518			hev_.vr = vr;
519
520			hydjetTree_->Fill();
521
522			}
523
524
525			// ------------ method called once each job just before starting event loop ------------
526			void
527			HeavyIonJetAnalyzer::beginJob(const edm::EventSetup& iSetup)
528			{
529	yilmaz	1.4	cout<<"Begin"<<endl;
530
531	yilmaz	1.1	iSetup.getData(pdt);
532
533			if(printLists_){
534			out_b.open(fBFileName.c_str());
535			if(out_b.good() == false)
536			throw cms::Exception("BadFile") << "Can\'t open file " << fBFileName;
537			out_n.open(fNFileName.c_str());
538			if(out_n.good() == false)
539			throw cms::Exception("BadFile") << "Can\'t open file " << fNFileName;
540			out_m.open(fMFileName.c_str());
541			if(out_m.good() == false)
542			throw cms::Exception("BadFile") << "Can\'t open file " << fMFileName;
543			}
544
545			if(doAnalysis_){
546			hydjetTree_ = f->make<TTree>("hi","Tree of Hydjet Events");
547	yilmaz	1.4	jetTree_ = f->make<TTree>("jet","Tree of Jet Constituents");
548	yilmaz	1.1
549			hydjetTree_->Branch("event",&hev_.event,"event/I");
550			hydjetTree_->Branch("b",&hev_.b,"b/F");
551			hydjetTree_->Branch("npart",&hev_.npart,"npart/F");
552			hydjetTree_->Branch("ncoll",&hev_.ncoll,"ncoll/F");
553			hydjetTree_->Branch("nhard",&hev_.nhard,"nhard/F");
554	yilmaz	1.4	hydjetTree_->Branch("phi0",&hev_.phi0,"phi0/F");
555	yilmaz	1.1
556			if(doParticles_){
557			hydjetTree_->Branch("n",hev_.n,"n[3]/I");
558			hydjetTree_->Branch("ptav",hev_.ptav,"ptav[3]/F");
559			hydjetTree_->Branch("np",&hev_.np,"np/I");
560			hydjetTree_->Branch("par_pt",hev_.par_pt,"par_pt[np]/F");
561			hydjetTree_->Branch("par_eta",hev_.par_eta,"par_eta[np]/F");
562			hydjetTree_->Branch("par_phi",hev_.par_phi,"par_phi[np]/F");
563			hydjetTree_->Branch("pdg",hev_.pdg,"pdg[np]/I");
564			hydjetTree_->Branch("chg",hev_.chg,"chg[np]/I");
565			}
566
567			hydjetTree_->Branch("vx",&hev_.vx,"vx/F");
568			hydjetTree_->Branch("vy",&hev_.vy,"vy/F");
569			hydjetTree_->Branch("vz",&hev_.vz,"vz/F");
570			hydjetTree_->Branch("vr",&hev_.vr,"vr/F");
571
572			hydjetTree_->Branch("njet",&hev_.njet,"njet/I");
573			hydjetTree_->Branch("et",hev_.et,"et[njet]/F");
574			hydjetTree_->Branch("eta",hev_.eta,"eta[njet]/F");
575			hydjetTree_->Branch("phi",hev_.phi,"phi[njet]/F");
576			hydjetTree_->Branch("area",hev_.area,"area[njet]/F");
577
578	yilmaz	1.4	hydjetTree_->Branch("r20",hev_.r20,"r20[njet]/F");
579			hydjetTree_->Branch("r50",hev_.r50,"r50[njet]/F");
580			hydjetTree_->Branch("r90",hev_.r90,"r90[njet]/F");
581
582			hydjetTree_->Branch("e01",hev_.e01,"e01[njet]/F");
583			hydjetTree_->Branch("e02",hev_.e02,"e02[njet]/F");
584			hydjetTree_->Branch("e03",hev_.e03,"e03[njet]/F");
585			hydjetTree_->Branch("e04",hev_.e04,"e04[njet]/F");
586			hydjetTree_->Branch("e05",hev_.e05,"e05[njet]/F");
587
588			hydjetTree_->Branch("ncons",hev_.ncons,"ncons[njet]/I");
589
590			jetTree_->Branch("event",&hev_.event,"event/I");
591			jetTree_->Branch("b",&hev_.b,"b/F");
592
593			jetTree_->Branch("ncons",&yet_.ncons,"ncons/I");
594			jetTree_->Branch("et",yet_.et,"et[ncons]/F");
595			jetTree_->Branch("eta",yet_.eta,"eta[ncons]/F");
596			jetTree_->Branch("phi",yet_.phi,"phi[ncons]/F");
597			jetTree_->Branch("dr",yet_.dr,"dr[ncons]/F");
598
599			jetTree_->Branch("etjet",&yet_.etjet,"etjet/F");
600			jetTree_->Branch("etajet",&yet_.etajet,"etajet/F");
601			jetTree_->Branch("phijet",&yet_.phijet,"phijet/F");
602
603			jetTree_->Branch("area",&yet_.area,"area/F");
604			jetTree_->Branch("r20",&yet_.r20,"r20/F");
605			jetTree_->Branch("r50",&yet_.r50,"r50/F");
606			jetTree_->Branch("r90",&yet_.r90,"r90/F");
607
608			jetTree_->Branch("e01",&yet_.e01,"e01/F");
609			jetTree_->Branch("e02",&yet_.e02,"e02/F");
610			jetTree_->Branch("e03",&yet_.e03,"e03/F");
611			jetTree_->Branch("e04",&yet_.e04,"e04/F");
612			jetTree_->Branch("e05",&yet_.e05,"e05/F");
613
614
615
616	yilmaz	1.1	}
617
618			}
619
620			// ------------ method called once each job just after ending the event loop ------------
621			void
622			HeavyIonJetAnalyzer::endJob() {
623			}
624
625			//define this as a plug-in
626			DEFINE_FWK_MODULE(HeavyIonJetAnalyzer);