ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/Jeng/PVStudy/plugins/PVStudy.cc

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.5 by jengbou, Wed Sep 9 05:56:50 2009 UTC vs.
Revision 1.6 by jengbou, Wed Sep 9 23:05:38 2009 UTC

#	Line 87 | Line 87 | PVStudy::PVStudy(const edm::ParameterSet
87
88		setRootStyle();
89
90	<	h["nTrk"]           = subDir.make<TH1F>("nTrk","Num of total tracks",300,0,300);
91	<	h["nTrkPV"]         = subDir.make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
92	<	h["trkPt"]          = subDir.make<TH1D>("trkPt","Pt of all tracks",100,0,100);
93	<	h["trkPtPV"]        = subDir.make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);
90	>	h["nTrk"]           = subDir.make<TH1F>("nTrk","Num of total rec tracks",300,0,300);
91	>	h["nTrkPV"]         = subDir.make<TH1F>("nTrkPV","Num of rec tracks of PVs",300,0,300);
92	>	h["trkPt"]          = subDir.make<TH1D>("trkPt","Pt of all rec tracks",100,0,100);
93	>	h["trkPtPV"]        = subDir.make<TH1D>("trkPtPV","Pt dist of rec tracks from PV",100,0,100);
94	>	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
95	>	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
96	>	h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
97	>	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
98
99		if (!realData_) {
100		h["deltax"]       = subDir.make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
#	Line 210 | Line 214 | std::vector<PVStudy::simPrimaryVertex> P
214		cout << "number of vertices " << evt->vertices_size() << endl;
215		}
216
217	<	int idx=0;
217	>	int idx=0; int npv=0;
218		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
219	<	vitr != evt->vertices_end(); ++vitr )
220	<	{ // loop for vertex ...
221	<	HepMC::FourVector pos = (*vitr)->position();
222	<	//HepLorentzVector pos = (*vitr)->position();
223	<	if (pos.t()>0) { continue;} // ?
219	>	vitr != evt->vertices_end(); ++vitr ) { // loop for vertex ...
220	>	HepMC::FourVector pos = (*vitr)->position();
221	>	//HepLorentzVector pos = (*vitr)->position();
222	>
223	>	// t component of PV:
224	>	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
225	>	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
226	>	//        std::cout << "Status = " << (*mother)->status() << std::endl;
227	>	HepMC::GenVertex * mv=(*mother)->production_vertex();
228	>	if( ((*mother)->status() == 3) && (!mv)) {
229	>	//        std::cout << "npv= " << npv << std::endl;
230	>	if (npv == 0) {
231	>	h["genPart_cT"]->Fill(pos.t()); // mm
232	>	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
233	>	}
234	>	npv++;
235	>	}
236	>	}
237	>	//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
238
239	<	bool hasMotherVertex=false;
240	<	//      if(verbose_) cout << "mothers" << endl;
241	<	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
242	<	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
243	<	HepMC::GenVertex * mv=(*mother)->production_vertex();
244	<	if (mv) {hasMotherVertex=true;}
245	<	//      if(verbose_) {
246	<	//        cout << "\t";
247	<	//        (*mother)->print();
230	<	//      }
239	>	bool hasMotherVertex=false;
240	>	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
241	>	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
242	>	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
243	>	HepMC::GenVertex * mv=(*mother)->production_vertex();
244	>	//        std::cout << "Status = " << (*mother)->status() << std::endl;
245	>	if (mv) {
246	>	hasMotherVertex=true;
247	>	if(!verbose_) break; //if verbose_, print all particles of gen vertices
248		}
249	+	if(verbose_) {
250	+	cout << "\t";
251	+	(*mother)->print();
252	+	}
253	+	}
254
255	<	if(hasMotherVertex) {continue;}
255	>	if(hasMotherVertex) continue;
256
257	<	// could be a new vertex, check  all primaries found so far to avoid multiple entries
258	<	const double mm2cm=0.1;
259	<	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
260	<	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
261	<	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
262	<	v0!=simpv.end(); v0++){
263	<	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
264	<	vp=&(*v0);
265	<	break;
244	<	}
257	>	// could be a new vertex, check  all primaries found so far to avoid multiple entries
258	>	const double mm2cm=0.1;
259	>	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
260	>	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
261	>	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
262	>	v0!=simpv.end(); v0++){
263	>	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
264	>	vp=&(*v0);
265	>	break;
266		}
267	+	}
268
269	<	if(!vp){
270	<	// this is a new vertex
271	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
272	<	simpv.push_back(sv);
273	<	vp=&simpv.back();
274	<	}else{
275	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
276	<	}
277	<	vp->genVertex.push_back((*vitr)->barcode());
278	<	// collect final state descendants
279	<	for ( HepMC::GenVertex::particle_iterator
280	<	daughter  = (*vitr)->particles_begin(HepMC::descendants);
281	<	daughter != (*vitr)->particles_end(HepMC::descendants);
282	<	++daughter ) {
283	<	if (isFinalstateParticle(*daughter)){
284	<	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
285	<	== vp->finalstateParticles.end()){
286	<	vp->finalstateParticles.push_back((*daughter)->barcode());
287	<	HepMC::FourVector m=(*daughter)->momentum();
288	<	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
289	<	// but adding FourVectors seems not to be foreseen
290	<	vp->ptot.setPx(vp->ptot.px()+m.px());
291	<	vp->ptot.setPy(vp->ptot.py()+m.py());
292	<	vp->ptot.setPz(vp->ptot.pz()+m.pz());
293	<	vp->ptot.setE(vp->ptot.e()+m.e());
294	<	vp->ptsq+=(m.perp())*(m.perp());
295	<	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
274	<	vp->nGenTrk++;
275	<	}
269	>	if(!vp){
270	>	// this is a new vertex
271	>	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
272	>	simpv.push_back(sv);
273	>	vp=&simpv.back();
274	>	}else{
275	>	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
276	>	}
277	>	vp->genVertex.push_back((*vitr)->barcode());
278	>	// collect final state descendants
279	>	for ( HepMC::GenVertex::particle_iterator daughter  = (*vitr)->particles_begin(HepMC::descendants);
280	>	daughter != (*vitr)->particles_end(HepMC::descendants);
281	>	++daughter ) {
282	>	if (isFinalstateParticle(*daughter)){
283	>	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
284	>	== vp->finalstateParticles.end()){
285	>	vp->finalstateParticles.push_back((*daughter)->barcode());
286	>	HepMC::FourVector m=(*daughter)->momentum();
287	>	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
288	>	// but adding FourVectors seems not to be foreseen
289	>	vp->ptot.setPx(vp->ptot.px()+m.px());
290	>	vp->ptot.setPy(vp->ptot.py()+m.py());
291	>	vp->ptot.setPz(vp->ptot.pz()+m.pz());
292	>	vp->ptot.setE(vp->ptot.e()+m.e());
293	>	vp->ptsq+=(m.perp())*(m.perp());
294	>	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
295	>	vp->nGenTrk++;
296		}
297		}
298		}
279	–	idx++;
299		}
300	+	}
301		}
302		return simpv;
303		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines