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.4 2009/06/16 12:34:34 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 "SimDataFormats/Vertex/interface/SimVertex.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.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_;

   vector<string> jetSrc_;

   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<vector<string> >("jetSrc");

   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("source",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_[0],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_){   
   edm::Handle<edm::SimVertexContainer> simVertices;
   iEvent.getByType<edm::SimVertexContainer>(simVertices);
   
   if (! simVertices.isValid() ) throw cms::Exception("FatalError") << "No vertices found\n";
   int inum = 0;
   
   edm::SimVertexContainer::const_iterator it=simVertices->begin();
   SimVertex vertex = (*it);
   cout<<" Vertex position "<< inum <<" " << vertex.position().rho()<<" "<<vertex.position().z()<<endl;
   vx = vertex.position().x();
   vy = vertex.position().y();
   vz = vertex.position().z();
   vr = vertex.position().rho();
   }

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