CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DoCastorAnalysis.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Forward
4 // Class : DoCastorAnalysis
5 //
6 // Implementation:
7 // <Notes on implementation>
8 //
9 // Original Author: P. Katsas
10 // Created: 02/2007
11 //
15 
18 
19 #include "TFile.h"
20 #include <cmath>
21 #include <iostream>
22 #include <iomanip>
23 #include <memory>
24 #include <vector>
25 
26 #define debug 0
27 
29 
30  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("DoCastorAnalysis");
31  verbosity = m_Anal.getParameter<int>("Verbosity");
32 
33  TreeFileName = m_Anal.getParameter<std::string>("CastorTreeFileName");
34 
35  if (verbosity > 0) {
36 
37  std::cout<<std::endl;
38  std::cout<<"============================================================================"<<std::endl;
39  std::cout << "DoCastorAnalysis:: Initialized as observer"<< std::endl;
40 
41  std::cout <<" Castor Tree will be created"<< std::endl;
42  std::cout <<" Castor Tree will be in file: "<<TreeFileName<<std::endl;
43  if(debug) getchar();
44 
45  std::cout<<"============================================================================"<<std::endl;
46  std::cout<<std::endl;
47  }
48 
49  std::cout << "DoCastorAnalysis: output event root file created"<< std::endl;
50  TString treefilename = TreeFileName;
51  CastorOutputEventFile = new TFile(treefilename,"RECREATE");
52 
53  CastorTree = new TTree("Sim","Sim");
54 
55  CastorTree->Branch("simhit_x","std::vector<double>",&psimhit_x);
56  CastorTree->Branch("simhit_y","std::vector<double>",&psimhit_y);
57  CastorTree->Branch("simhit_z","std::vector<double>",&psimhit_z);
58 
59  CastorTree->Branch("simhit_eta","std::vector<double>",&psimhit_eta);
60  CastorTree->Branch("simhit_phi","std::vector<double>",&psimhit_phi);
61  CastorTree->Branch("simhit_energy","std::vector<double>",&psimhit_energy);
62 
63  // CastorTree->Branch("simhit_time","std::vector<double>",&psimhit_time);
64  CastorTree->Branch("simhit_sector","std::vector<int>",&psimhit_sector);
65  CastorTree->Branch("simhit_module","std::vector<int>",&psimhit_module);
66 
67  CastorTree->Branch("simhit_etot",&simhit_etot,"simhit_etot/D");
68 }
69 
71 
72  //destructor of DoCastorAnalysis
73 
75  //-- CastorOutputEventFile->Write();
76  CastorTree->Write("",TObject::kOverwrite);
77  std::cout << "DoCastorAnalysis: Ntuple event written" << std::endl;
78  if(debug) getchar();
79  CastorOutputEventFile->Close();
80  std::cout << "DoCastorAnalysis: Event file closed" << std::endl;
81  if(debug) getchar();
82 
83  if (verbosity > 0) {
84  std::cout<<std::endl<<"DoCastorAnalysis: end of process"<<std::endl;
85  if(debug) getchar();
86  }
87 
88 }
89 
90 //=================================================================== per EVENT
91 
93 
94  std::cout << " Starting new job " << std::endl;
95 }
96 
97 //==================================================================== per RUN
98 
100 
101  std::cout << std::endl << "DoCastorAnalysis: Starting Run"<< std::endl;
102 
103  // std::cout << "DoCastorAnalysis: output event root file created"<< std::endl;
104  // TString treefilename = TreeFileName;
105  // CastorOutputEventFile = new TFile(treefilename,"RECREATE");
106 
107  eventIndex = 1;
108 }
109 
111  std::cout << "DoCastorAnalysis: Processing Event Number: "<<eventIndex<< std::endl;
112  eventIndex++;
113 }
114 
115 
116 //================= End of EVENT ===============
117 
119 
120  // Look for the Hit Collection
121 
122  // access to the G4 hit collections
123  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
124 
125  int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("CastorFI");
126  CaloG4HitCollection* theCAFI = (CaloG4HitCollection*) allHC->GetHC(CAFIid);
127 
128  CastorNumberingScheme *theCastorNumScheme = new CastorNumberingScheme();
129 
130  unsigned int volumeID=0;
131  // std::map<int,float,std::less<int> > themap;
132 
133  int nentries = theCAFI->entries();
134  if(debug) std::cout<<"nentries in CAFI: "<<nentries<<std::endl;
135  if(debug) getchar();
136 
138  psimhit_x->clear();
139  psimhit_x->reserve(nentries);
140 
142  psimhit_y->clear();
143  psimhit_y->reserve(nentries);
144 
146  psimhit_z->clear();
147  psimhit_z->reserve(nentries);
148 
150  psimhit_eta->clear();
151  psimhit_eta->reserve(nentries);
152 
154  psimhit_phi->clear();
155  psimhit_phi->reserve(nentries);
156 
158  psimhit_energy->clear();
159  psimhit_energy->reserve(nentries);
160 
161  //psimhit_time=&simhit_time;
162  //psimhit_time->clear();
163  //psimhit_time->reserve(nentries);
164 
166  psimhit_sector->clear();
167  psimhit_sector->reserve(nentries);
168 
170  psimhit_module->clear();
171  psimhit_module->reserve(nentries);
172 
173  simhit_etot = 0;
174 
175  if (nentries > 0) {
176 
177  for (int ihit = 0; ihit < nentries; ihit++) {
178  CaloG4Hit* aHit = (*theCAFI)[ihit];
179  volumeID = aHit->getUnitID();
180 
181  //themap[volumeID] += aHit->getEnergyDeposit();
182  int zside,sector,zmodule;
183 
184  theCastorNumScheme->unpackIndex(volumeID,zside,sector,zmodule);
185 
186  double energy = aHit->getEnergyDeposit()/GeV;
187  //double time = aHit->getTimeSlice();
188 
189  math::XYZPoint pos = aHit->getPosition();
190  double theta = pos.theta();
191  double eta = -log(tan(theta/2.));
192  double phi = pos.phi();
193 
194  psimhit_x->push_back(pos.x());
195  psimhit_y->push_back(pos.y());
196  psimhit_z->push_back(pos.z());
197 
198  psimhit_eta->push_back(eta);
199  psimhit_phi->push_back(phi);
200  psimhit_energy->push_back(energy);
201 
202  // psimhit_time->push_back(time);
203  psimhit_sector->push_back(sector);
204  psimhit_module->push_back(zmodule);
205 
207 
208  if(debug) std::cout<<"hit "<<ihit+1<<" : x = "<<(*psimhit_x)[ihit]<<" , eta = "<<(*psimhit_eta)[ihit]
209  <<" , phi = "<<(*psimhit_phi)[ihit]<<" , energy = "<<(*psimhit_energy)[ihit]<<std::endl;
210  }
211 
212  //if(debug) std::cout<<" total energy = "<<simhit_etot<<std::endl;
213  if(debug) getchar();
214  CastorTree->Fill();
215 
216  } // nentries > 0
217 }
218 
220 
221 void DoCastorAnalysis::update(const G4Step * aStep) {;}
222 
223 
std::vector< double > * psimhit_z
T getParameter(std::string const &) const
math::XYZPoint getPosition() const
Definition: CaloG4Hit.h:56
std::vector< double > simhit_eta
std::vector< double > * psimhit_x
std::vector< double > simhit_x
std::vector< double > * psimhit_energy
std::vector< double > * psimhit_phi
Geom::Theta< T > theta() const
DoCastorAnalysis(const edm::ParameterSet &p)
virtual ~DoCastorAnalysis()
T eta() const
std::vector< double > simhit_phi
static void unpackIndex(const uint32_t &idx, int &z, int &sector, int &zmodule)
std::vector< double > * psimhit_y
std::vector< double > simhit_y
std::vector< double > * psimhit_eta
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
TFile * CastorOutputEventFile
void update(const BeginOfJob *run)
This routine will be called when the appropriate signal arrives.
std::vector< int > simhit_sector
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
std::vector< int > * psimhit_sector
std::string TreeFileName
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
tuple job
Definition: launcher.py:57
std::vector< double > simhit_energy
tuple cout
Definition: gather_cfg.py:121
uint32_t getUnitID() const
Definition: CaloG4Hit.h:69
#define debug
Definition: MEtoEDMFormat.h:34
std::vector< int > * psimhit_module
double getEnergyDeposit() const
Definition: CaloG4Hit.h:81
std::vector< int > simhit_module
std::vector< double > simhit_z
Definition: DDAxes.h:10