hpstr
The Heavy Photon Search Toolkit for Reconstruction (hpstr) provides an interface to physics data from the HPS experiment saved in the LCIO format and converts it into an ROOT based format.
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
VertexAnaProcessor.cxx
Go to the documentation of this file.
1
7#include "VertexAnaProcessor.h"
8#include <iostream>
9#include <fstream>
10#include <map>
11
12VertexAnaProcessor::VertexAnaProcessor(const std::string& name, Process& process) : Processor(name,process) {
13
14}
15
16//TODO Check this destructor
17
18VertexAnaProcessor::~VertexAnaProcessor(){}
19
20void VertexAnaProcessor::configure(const ParameterSet& parameters) {
21 std::cout << "Configuring VertexAnaProcessor" <<std::endl;
22 try
23 {
24 debug_ = parameters.getInteger("debug",debug_);
25 anaName_ = parameters.getString("anaName",anaName_);
26 tsColl_ = parameters.getString("tsColl",tsColl_);
27 vtxColl_ = parameters.getString("vtxColl",vtxColl_);
28 trkColl_ = parameters.getString("trkColl",trkColl_);
29 hitColl_ = parameters.getString("hitColl",hitColl_);
30 ecalColl_ = parameters.getString("ecalColl",ecalColl_);
31 mcColl_ = parameters.getString("mcColl",mcColl_);
32 isRadPDG_ = parameters.getInteger("isRadPDG",isRadPDG_);
33 makeFlatTuple_ = parameters.getInteger("makeFlatTuple",makeFlatTuple_);
34
35 selectionCfg_ = parameters.getString("vtxSelectionjson",selectionCfg_);
36 histoCfg_ = parameters.getString("histoCfg",histoCfg_);
37 mcHistoCfg_ = parameters.getString("mcHistoCfg",mcHistoCfg_);
38 timeOffset_ = parameters.getDouble("CalTimeOffset",timeOffset_);
39 beamE_ = parameters.getDouble("beamE",beamE_);
40 isData_ = parameters.getInteger("isData",isData_);
41 analysis_ = parameters.getString("analysis");
42
43 //region definitions
44 regionSelections_ = parameters.getVString("regionDefinitions",regionSelections_);
45
46 //v0 projection fits
47 v0ProjectionFitsCfg_ = parameters.getString("v0ProjectionFitsCfg", v0ProjectionFitsCfg_);
48
49 //beamspot positions
50 beamPosCfg_ = parameters.getString("beamPosCfg", beamPosCfg_);
51 //track time bias corrections
52 eleTrackTimeBias_ = parameters.getDouble("eleTrackTimeBias",eleTrackTimeBias_);
53 posTrackTimeBias_ = parameters.getDouble("posTrackTimeBias",posTrackTimeBias_);
54
55 }
56 catch (std::runtime_error& error)
57 {
58 std::cout<<error.what()<<std::endl;
59 }
60}
61
62void VertexAnaProcessor::initialize(TTree* tree) {
63 tree_ = tree;
64 _ah = std::make_shared<AnaHelpers>();
65
66 vtxSelector = std::make_shared<BaseSelector>(anaName_+"_"+"vtxSelection",selectionCfg_);
67 vtxSelector->setDebug(debug_);
68 vtxSelector->LoadSelection();
69
70 _vtx_histos = std::make_shared<TrackHistos>(anaName_+"_"+"vtxSelection");
71 _vtx_histos->loadHistoConfig(histoCfg_);
72 _vtx_histos->DefineHistos();
73
74 if(!isData_ && mc_reg_on_){
75 _mc_vtx_histos = std::make_shared<MCAnaHistos>(anaName_+"_mc_"+"vtxSelection");
76 _mc_vtx_histos->loadHistoConfig(mcHistoCfg_);
77 _mc_vtx_histos->DefineHistos();
78 _mc_vtx_histos->Define2DHistos();
79 }
80
81 //Load Run Dependent V0 target projection fits from json
82 if(!v0ProjectionFitsCfg_.empty()){
83 std::ifstream v0proj_file(v0ProjectionFitsCfg_);
84 v0proj_file >> v0proj_fits_;
85 v0proj_file.close();
86 }
87
88 //Run Dependent Corrections
89 //Beam Position
90 if(!beamPosCfg_.empty()){
91 std::ifstream bpc_file(beamPosCfg_);
92 bpc_file >> bpc_configs_;
93 bpc_file.close();
94 }
95 // histos = new MCAnaHistos(anaName_);
96 //histos->loadHistoConfig(histCfgFilename_)
97 //histos->DefineHistos();
98 //histos->Define2DHistos();
99
100
101 //For each region initialize plots
102
103 for (unsigned int i_reg = 0; i_reg < regionSelections_.size(); i_reg++) {
104 std::string regname = AnaHelpers::getFileName(regionSelections_[i_reg],false);
105 std::cout<<"Setting up region:: " << regname <<std::endl;
106 _reg_vtx_selectors[regname] = std::make_shared<BaseSelector>(anaName_+"_"+regname, regionSelections_[i_reg]);
107 _reg_vtx_selectors[regname]->setDebug(debug_);
108 _reg_vtx_selectors[regname]->LoadSelection();
109
110 _reg_vtx_histos[regname] = std::make_shared<TrackHistos>(anaName_+"_"+regname);
111 _reg_vtx_histos[regname]->loadHistoConfig(histoCfg_);
112 _reg_vtx_histos[regname]->DefineHistos();
113
114
115 bool mc_reg_on_ = false;
116 if(!isData_ && mc_reg_on_){
117 _reg_mc_vtx_histos[regname] = std::make_shared<MCAnaHistos>(anaName_+"_mc_"+regname);
118 _reg_mc_vtx_histos[regname]->loadHistoConfig(mcHistoCfg_);
119 _reg_mc_vtx_histos[regname]->DefineHistos();
120 }
121
122 //Build a flat tuple for vertex and track params
123 if (makeFlatTuple_){
124 _reg_tuples[regname] = std::make_shared<FlatTupleMaker>(anaName_+"_"+regname+"_tree");
125
126 //vtx vars
127 _reg_tuples[regname]->addVariable("unc_vtx_mass");
128 _reg_tuples[regname]->addVariable("unc_vtx_z");
129 _reg_tuples[regname]->addVariable("unc_vtx_chi2");
130 _reg_tuples[regname]->addVariable("unc_vtx_psum");
131 _reg_tuples[regname]->addVariable("unc_vtx_px");
132 _reg_tuples[regname]->addVariable("unc_vtx_py");
133 _reg_tuples[regname]->addVariable("unc_vtx_pz");
134 _reg_tuples[regname]->addVariable("unc_vtx_x");
135 _reg_tuples[regname]->addVariable("unc_vtx_y");
136 _reg_tuples[regname]->addVariable("unc_vtx_ele_pos_clus_dt");
137 _reg_tuples[regname]->addVariable("run_number");
138 _reg_tuples[regname]->addVariable("unc_vtx_cxx");
139 _reg_tuples[regname]->addVariable("unc_vtx_cyy");
140 _reg_tuples[regname]->addVariable("unc_vtx_czz");
141 _reg_tuples[regname]->addVariable("unc_vtx_cyx");
142 _reg_tuples[regname]->addVariable("unc_vtx_czy");
143 _reg_tuples[regname]->addVariable("unc_vtx_czx");
144 _reg_tuples[regname]->addVariable("unc_vtx_proj_x");
145 _reg_tuples[regname]->addVariable("unc_vtx_proj_y");
146 _reg_tuples[regname]->addVariable("unc_vtx_proj_x_sig");
147 _reg_tuples[regname]->addVariable("unc_vtx_proj_y_sig");
148 _reg_tuples[regname]->addVariable("unc_vtx_proj_sig");
149
150 //track vars
151 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_p");
152 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_t");
153 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_d0");
154 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_phi0");
155 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_omega");
156 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_tanLambda");
157 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_z0");
158 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_chi2ndf");
159 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_clust_dt");
160 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_z0Err");
161 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_d0Err");
162 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_tanLambdaErr");
163 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_PhiErr");
164 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_OmegaErr");
165 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_L1_isolation");
166 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_nhits");
167 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_x");
168 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_y");
169 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_z");
170 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_px");
171 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_py");
172 _reg_tuples[regname]->addVariable("unc_vtx_ele_track_pz");
173
174 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_clust_dt");
175 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_p");
176 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_t");
177 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_d0");
178 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_phi0");
179 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_omega");
180 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_tanLambda");
181 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_z0");
182 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_chi2ndf");
183 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_z0Err");
184 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_d0Err");
185 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_tanLambdaErr");
186 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_PhiErr");
187 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_OmegaErr");
188 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_L1_isolation");
189 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_nhits");
190 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_x");
191 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_y");
192 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_z");
193 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_px");
194 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_py");
195 _reg_tuples[regname]->addVariable("unc_vtx_pos_track_pz");
196
197 //clust vars
198 _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_E");
199 _reg_tuples[regname]->addVariable("unc_vtx_ele_clust_corr_t");
200
201 _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_E");
202 _reg_tuples[regname]->addVariable("unc_vtx_pos_clust_corr_t");
203
204 if(!isData_)
205 {
206 _reg_tuples[regname]->addVariable("true_vtx_z");
207 _reg_tuples[regname]->addVariable("true_vtx_mass");
208 _reg_tuples[regname]->addVariable("ap_true_vtx_z");
209 _reg_tuples[regname]->addVariable("ap_true_vtx_mass");
210 _reg_tuples[regname]->addVariable("ap_true_vtx_energy");
211 _reg_tuples[regname]->addVariable("vd_true_vtx_z");
212 _reg_tuples[regname]->addVariable("vd_true_vtx_mass");
213 _reg_tuples[regname]->addVariable("vd_true_vtx_energy");
214 _reg_tuples[regname]->addVariable("hitCode");
215 _reg_tuples[regname]->addVariable("L1hitCode");
216 _reg_tuples[regname]->addVariable("L2hitCode");
217 }
218 }
219
220 _regions.push_back(regname);
221 }
222
223 // Get list of branches in tree to help protect accessing them
224 int nBr = tree_->GetListOfBranches()->GetEntries();
225 if (debug_) std::cout << "Tree has " << nBr << " branches" << std::endl;
226 for(int iBr = 0; iBr < nBr; iBr++)
227 {
228 TBranch *br = dynamic_cast<TBranch*>(tree_->GetListOfBranches()->At(iBr));
229 brMap_.insert(std::map<const char *, int, char_cmp>::value_type(br->GetName(), 1));
230 if (debug_) std::cout << br->GetName() << ": " << brMap_[br->GetName()] << std::endl;
231 }
232
233 //init Reading Tree
234 tree_->SetBranchAddress("EventHeader", &evth_ , &bevth_);
235 if (brMap_.find(tsColl_.c_str()) != brMap_.end()) tree_->SetBranchAddress(tsColl_.c_str(), &ts_ , &bts_);
236 tree_->SetBranchAddress(vtxColl_.c_str(), &vtxs_ , &bvtxs_);
237 if (brMap_.find(hitColl_.c_str()) != brMap_.end()) tree_->SetBranchAddress(hitColl_.c_str(), &hits_ , &bhits_);
238 tree_->SetBranchAddress(ecalColl_.c_str(), &ecal_ , &becal_);
239 if(!isData_ && !mcColl_.empty()) tree_->SetBranchAddress(mcColl_.c_str() , &mcParts_, &bmcParts_);
240 //If track collection name is empty take the tracks from the particles. TODO:: change this
241 if (!trkColl_.empty())
242 tree_->SetBranchAddress(trkColl_.c_str(),&trks_, &btrks_);
243}
244
245bool VertexAnaProcessor::process(IEvent* ievent) {
246 if(debug_) {
247 std:: cout << "----------------- Event " << evth_->getEventNumber() << " -----------------" << std::endl;
248 }
249 HpsEvent* hps_evt = (HpsEvent*) ievent;
250 double weight = 1.;
251 int run_number = evth_->getRunNumber();
252 int closest_run;
253 if(!bpc_configs_.empty()){
254 for(auto run : bpc_configs_.items()){
255 int check_run = std::stoi(run.key());
256 if(check_run > run_number)
257 break;
258 else{
259 closest_run = check_run;
260 }
261 }
262 beamPosCorrections_ = {bpc_configs_[std::to_string(closest_run)]["beamspot_x"],
263 bpc_configs_[std::to_string(closest_run)]["beamspot_y"],
264 bpc_configs_[std::to_string(closest_run)]["beamspot_z"]};
265 }
266
267
268 //Get "true" values
269 //AP
270 double apMass = -0.9;
271 double apZ = -0.9;
272 double apEnergy = -0.9;
273 //Simp
274 double vdMass = -0.9;
275 double vdZ = -0.9;
276 double vdEnergy = -0.9;
277
278 //Plot info about which trigger bits are present in the event
279 if (ts_ != nullptr)
280 {
281 _vtx_histos->Fill2DHisto("trig_count_hh",
282 ((int)ts_->prescaled.Single_3_Top)+((int)ts_->prescaled.Single_3_Bot),
283 ((int)ts_->prescaled.Single_2_Top)+((int)ts_->prescaled.Single_2_Bot));
284 }
285 int NposTrks = 0;
286 int NeleTrks = 0;
287 for (int iT = 0; iT < trks_->size(); iT++)
288 {
289 if (trks_->at(iT)->getCharge() > 0) NposTrks++;
290 else NeleTrks++;
291 }
292 _vtx_histos->Fill2DHisto("n_tracks_hh", NeleTrks, NposTrks);
293 _vtx_histos->Fill1DHisto("n_vtx_h", vtxs_->size());
294
295 if (mcParts_) {
296 for(int i = 0; i < mcParts_->size(); i++)
297 {
298 if(mcParts_->at(i)->getPDG() == 622)
299 {
300 apMass = mcParts_->at(i)->getMass();
301 apZ = mcParts_->at(i)->getVertexPosition().at(2);
302 apEnergy = mcParts_->at(i)->getEnergy();
303 }
304 if(mcParts_->at(i)->getPDG() == 625)
305 {
306 vdMass = mcParts_->at(i)->getMass();
307 vdZ = mcParts_->at(i)->getVertexPosition().at(2);
308 vdEnergy = mcParts_->at(i)->getEnergy();
309 }
310 }
311
312 if (!isData_ && mc_reg_on_) _mc_vtx_histos->FillMCParticles(mcParts_, analysis_);
313 }
314 //Store processed number of events
315 std::vector<Vertex*> selected_vtxs;
316 bool passVtxPresel = false;
317
318 // Fill some diagnostic histos
319 for ( int i_ecal = 0; i_ecal < ecal_->size(); i_ecal++ ) {
320
321 if (vtxs_->size() == 0){
322 _vtx_histos->Fill1DHisto("EecalClus_noVtxs_h",ecal_->at(i_ecal)->getEnergy());
323 } else {
324 _vtx_histos->Fill1DHisto("EecalClus_isVtxs_h",ecal_->at(i_ecal)->getEnergy());
325 }
326 }
327
328
329 if (vtxs_->size() == 0){
330 _vtx_histos->Fill1DHisto("n_ecalClus_noVtxs_h",ecal_->size());
331 _vtx_histos->Fill1DHisto("n_tracks_noVtxs_h",trks_->size());
332 for (int i_trk = 0; i_trk < trks_->size(); i_trk++ ){
333 _vtx_histos->Fill1DHisto("Ptracks_noVtxs_h",trks_->at(i_trk)->getP());
334 }
335
336 } else {
337 _vtx_histos->Fill1DHisto("n_ecalClus_isVtxs_h",ecal_->size());
338 _vtx_histos->Fill1DHisto("n_tracks_isVtxs_h",trks_->size());
339 for (int i_trk = 0; i_trk < trks_->size(); i_trk++ ){
340 _vtx_histos->Fill1DHisto("Ptracks_isVtxs_h",trks_->at(i_trk)->getP());
341 }
342 }
343
344 if(debug_){
345 std::cout<<"Number of vertices found in event: "<< vtxs_->size()<<std::endl;
346 }
347
348 // Loop over vertices in event and make selections
349 for ( int i_vtx = 0; i_vtx < vtxs_->size(); i_vtx++ ) {
350 vtxSelector->getCutFlowHisto()->Fill(0.,weight);
351
352 Vertex* vtx = vtxs_->at(i_vtx);
353 Particle* ele = nullptr;
354 Track* ele_trk = nullptr;
355 Particle* pos = nullptr;
356 Track* pos_trk = nullptr;
357
358 //Trigger requirement - *really hate* having to do it here for each vertex.
359
360 if (isData_) {
361 if (!vtxSelector->passCutEq("Pair1_eq",(int)evth_->isPair1Trigger(),weight))
362 break;
363 }
364
365 bool foundParts = _ah->GetParticlesFromVtx(vtx,ele,pos);
366 if (!foundParts) {
367 if(debug_) std::cout<<"VertexAnaProcessor::WARNING::Found vtx without ele/pos. Skip."<<std::endl;
368 continue;
369 }
370
371 if (!trkColl_.empty()) {
372 bool foundTracks = _ah->MatchToGBLTracks((ele->getTrack()).getID(),(pos->getTrack()).getID(),
373 ele_trk, pos_trk, *trks_);
374
375 if (!foundTracks) {
376 if(debug_) std::cout<<"VertexAnaProcessor::ERROR couldn't find ele/pos in the GBLTracks collection"<<std::endl;
377 continue;
378 }
379 }
380 else {
381 ele_trk = (Track*)ele->getTrack().Clone();
382 pos_trk = (Track*)pos->getTrack().Clone();
383 }
384
385 //Beam Position Corrections
386 ele_trk->applyCorrection("z0", beamPosCorrections_.at(1));
387 pos_trk->applyCorrection("z0", beamPosCorrections_.at(1));
388 //Track Time Corrections
389 ele_trk->applyCorrection("track_time",eleTrackTimeBias_);
390 pos_trk->applyCorrection("track_time", posTrackTimeBias_);
391
392 //Add the momenta to the tracks - do not do that
393 //ele_trk->setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]);
394 //pos_trk->setMomentum(pos->getMomentum()[0],pos->getMomentum()[1],pos->getMomentum()[2]);
395
396 double ele_E = ele->getEnergy();
397 double pos_E = pos->getEnergy();
398
399 CalCluster eleClus = ele->getCluster();
400 CalCluster posClus = pos->getCluster();
401
402
403 //Compute analysis variables here.
404 TLorentzVector p_ele;
405 p_ele.SetPxPyPzE(ele_trk->getMomentum()[0],ele_trk->getMomentum()[1],ele_trk->getMomentum()[2],ele->getEnergy());
406 TLorentzVector p_pos;
407 p_pos.SetPxPyPzE(pos_trk->getMomentum()[0],pos_trk->getMomentum()[1],pos_trk->getMomentum()[2],ele->getEnergy());
408
409 //Tracks in opposite volumes - useless
410 //if (!vtxSelector->passCutLt("eleposTanLambaProd_lt",ele_trk->getTanLambda() * pos_trk->getTanLambda(),weight))
411 // continue;
412
413 //Ele Track Time
414 if (!vtxSelector->passCutLt("eleTrkTime_lt",fabs(ele_trk->getTrackTime()),weight))
415 continue;
416
417 //Pos Track Time
418 if (!vtxSelector->passCutLt("posTrkTime_lt",fabs(pos_trk->getTrackTime()),weight))
419 continue;
420
421 //Ele Track-cluster match
422 if (!vtxSelector->passCutLt("eleTrkCluMatch_lt",ele->getGoodnessOfPID(),weight))
423 continue;
424
425 //Pos Track-cluster match
426 if (!vtxSelector->passCutLt("posTrkCluMatch_lt",pos->getGoodnessOfPID(),weight))
427 continue;
428
429 //Require Positron Cluster exists
430 if (!vtxSelector->passCutGt("posClusE_gt",posClus.getEnergy(),weight))
431 continue;
432
433 //Require Positron Cluster does NOT exists
434 if (!vtxSelector->passCutLt("posClusE_lt",posClus.getEnergy(),weight))
435 continue;
436
437
438 double corr_eleClusterTime = ele->getCluster().getTime() - timeOffset_;
439 double corr_posClusterTime = pos->getCluster().getTime() - timeOffset_;
440
441 double botClusTime = 0.0;
442 if(ele->getCluster().getPosition().at(1) < 0.0) botClusTime = ele->getCluster().getTime();
443 else botClusTime = pos->getCluster().getTime();
444
445 //Bottom Cluster Time
446 if (!vtxSelector->passCutLt("botCluTime_lt", botClusTime, weight))
447 continue;
448
449 if (!vtxSelector->passCutGt("botCluTime_gt", botClusTime, weight))
450 continue;
451
452 //Ele Pos Cluster Time Difference
453 if (!vtxSelector->passCutLt("eleposCluTimeDiff_lt",fabs(corr_eleClusterTime - corr_posClusterTime),weight))
454 continue;
455
456 //Ele Track-Cluster Time Difference
457 if (!vtxSelector->passCutLt("eleTrkCluTimeDiff_lt",fabs(ele_trk->getTrackTime() - corr_eleClusterTime),weight))
458 continue;
459
460 //Pos Track-Cluster Time Difference
461 if (!vtxSelector->passCutLt("posTrkCluTimeDiff_lt",fabs(pos_trk->getTrackTime() - corr_posClusterTime),weight))
462 continue;
463
464 TVector3 ele_mom;
465 //ele_mom.SetX(ele->getMomentum()[0]);
466 //ele_mom.SetY(ele->getMomentum()[1]);
467 //ele_mom.SetZ(ele->getMomentum()[2]);
468 ele_mom.SetX(ele_trk->getMomentum()[0]);
469 ele_mom.SetY(ele_trk->getMomentum()[1]);
470 ele_mom.SetZ(ele_trk->getMomentum()[2]);
471
472
473 TVector3 pos_mom;
474 //pos_mom.SetX(pos->getMomentum()[0]);
475 //pos_mom.SetY(pos->getMomentum()[1]);
476 //pos_mom.SetZ(pos->getMomentum()[2]);
477 pos_mom.SetX(pos_trk->getMomentum()[0]);
478 pos_mom.SetY(pos_trk->getMomentum()[1]);
479 pos_mom.SetZ(pos_trk->getMomentum()[2]);
480
481
482
483 //Ele Track Quality - Chi2
484 if (!vtxSelector->passCutLt("eleTrkChi2_lt",ele_trk->getChi2(),weight))
485 continue;
486
487 //Pos Track Quality - Chi2
488 if (!vtxSelector->passCutLt("posTrkChi2_lt",pos_trk->getChi2(),weight))
489 continue;
490
491 //Ele Track Quality - Chi2Ndf
492 if (!vtxSelector->passCutLt("eleTrkChi2Ndf_lt",ele_trk->getChi2Ndf(),weight))
493 continue;
494
495 //Pos Track Quality - Chi2Ndf
496 if (!vtxSelector->passCutLt("posTrkChi2Ndf_lt",pos_trk->getChi2Ndf(),weight))
497 continue;
498
499 //Beam Electron cut
500 if (!vtxSelector->passCutLt("eleMom_lt",ele_mom.Mag(),weight))
501 continue;
502
503 //Ele min momentum cut
504 if (!vtxSelector->passCutGt("eleMom_gt",ele_mom.Mag(),weight))
505 continue;
506
507 //Pos min momentum cut
508 if (!vtxSelector->passCutGt("posMom_gt",pos_mom.Mag(),weight))
509 continue;
510
511 //Ele nHits
512 int ele2dHits = ele_trk->getTrackerHitCount();
513 if (!ele_trk->isKalmanTrack())
514 ele2dHits*=2;
515
516 if (!vtxSelector->passCutGt("eleN2Dhits_gt",ele2dHits,weight)) {
517 continue;
518 }
519
520 //Pos nHits
521 int pos2dHits = pos_trk->getTrackerHitCount();
522 if (!pos_trk->isKalmanTrack())
523 pos2dHits*=2;
524
525 if (!vtxSelector->passCutGt("posN2Dhits_gt",pos2dHits,weight)) {
526 continue;
527 }
528
529 //Less than 4 shared hits for ele/pos track
530 if (!vtxSelector->passCutLt("eleNshared_lt",ele_trk->getNShared(),weight)) {
531 continue;
532 }
533
534 if (!vtxSelector->passCutLt("posNshared_lt",pos_trk->getNShared(),weight)) {
535 continue;
536 }
537
538
539 //Vertex Quality
540 if (!vtxSelector->passCutLt("chi2unc_lt",vtx->getChi2(),weight))
541 continue;
542
543 //Max vtx momentum
544 if (!vtxSelector->passCutLt("maxVtxMom_lt",(ele_mom+pos_mom).Mag(),weight))
545 continue;
546
547 //Min vtx momentum
548
549 if (!vtxSelector->passCutGt("minVtxMom_gt",(ele_mom+pos_mom).Mag(),weight))
550 continue;
551
552 _vtx_histos->Fill1DVertex(vtx,
553 ele,
554 pos,
555 ele_trk,
556 pos_trk,
557 weight);
558
559 double ele_pos_dt = corr_eleClusterTime - corr_posClusterTime;
560 double psum = ele_mom.Mag()+pos_mom.Mag();
561
562 _vtx_histos->Fill1DTrack(ele_trk,weight, "ele_");
563 _vtx_histos->Fill1DTrack(pos_trk,weight, "pos_");
564 _vtx_histos->Fill1DHisto("ele_track_n2dhits_h", ele2dHits, weight);
565 _vtx_histos->Fill1DHisto("pos_track_n2dhits_h", pos2dHits, weight);
566 _vtx_histos->Fill1DHisto("vtx_Psum_h", p_ele.P()+p_pos.P(), weight);
567 _vtx_histos->Fill1DHisto("vtx_Esum_h", ele_E + pos_E, weight);
568 _vtx_histos->Fill1DHisto("ele_pos_clusTimeDiff_h", (corr_eleClusterTime - corr_posClusterTime), weight);
569 _vtx_histos->Fill2DHisto("ele_vtxZ_iso_hh", TMath::Min(ele_trk->getIsolation(0), ele_trk->getIsolation(1)), vtx->getZ(), weight);
570 _vtx_histos->Fill2DHisto("pos_vtxZ_iso_hh", TMath::Min(pos_trk->getIsolation(0), pos_trk->getIsolation(1)), vtx->getZ(), weight);
571 _vtx_histos->Fill2DHistograms(vtx,weight);
572 _vtx_histos->Fill2DTrack(ele_trk,weight,"ele_");
573 _vtx_histos->Fill2DTrack(pos_trk,weight,"pos_");
574 _vtx_histos->Fill1DHisto("mcMass622_h",apMass);
575 _vtx_histos->Fill1DHisto("mcZ622_h",apZ);
576
577 //New SIMP histos for developing loose preselection cuts
578 //2d histos
579 _vtx_histos->Fill2DHisto("ele_clusT_v_ele_trackT_hh", ele_trk->getTrackTime(), corr_eleClusterTime, weight);
580 _vtx_histos->Fill2DHisto("pos_clusT_v_pos_trackT_hh", pos_trk->getTrackTime(), corr_posClusterTime, weight);
581 _vtx_histos->Fill2DHisto("ele_track_time_v_P_hh", ele_trk->getP(), ele_trk->getTrackTime(), weight);
582 _vtx_histos->Fill2DHisto("pos_track_time_v_P_hh", pos_trk->getP(), pos_trk->getTrackTime(), weight);
583 _vtx_histos->Fill2DHisto("ele_pos_clusTimeDiff_v_pSum_hh",ele_mom.Mag()+pos_mom.Mag(), ele_pos_dt, weight);
584 _vtx_histos->Fill2DHisto("ele_cluster_energy_v_track_p_hh",ele_trk->getP(), eleClus.getEnergy(), weight);
585 _vtx_histos->Fill2DHisto("pos_cluster_energy_v_track_p_hh",pos_trk->getP(), posClus.getEnergy(), weight);
586 _vtx_histos->Fill2DHisto("ele_track_cluster_dt_v_EoverP_hh",eleClus.getEnergy()/ele_trk->getP(), ele_trk->getTrackTime() - corr_eleClusterTime, weight);
587 _vtx_histos->Fill2DHisto("pos_track_cluster_dt_v_EoverP_hh",posClus.getEnergy()/pos_trk->getP(), pos_trk->getTrackTime() - corr_posClusterTime, weight);
588 _vtx_histos->Fill2DHisto("ele_track_clus_dt_v_p_hh",ele_trk->getP(), ele_trk->getTrackTime() - corr_eleClusterTime, weight);
589 _vtx_histos->Fill2DHisto("pos_track_clus_dt_v_p_hh",pos_trk->getP(), pos_trk->getTrackTime() - corr_posClusterTime, weight);
590 //chi2 2d plots
591 _vtx_histos->Fill2DHisto("ele_track_chi2ndf_v_time_hh", ele_trk->getTrackTime(), ele_trk->getChi2Ndf(), weight);
592 _vtx_histos->Fill2DHisto("ele_track_chi2ndf_v_p_hh", ele_trk->getP(), ele_trk->getChi2Ndf(), weight);
593 _vtx_histos->Fill2DHisto("ele_track_chi2ndf_v_tanlambda_hh", ele_trk->getTanLambda(), ele_trk->getChi2Ndf(), weight);
594 _vtx_histos->Fill2DHisto("ele_track_chi2ndf_v_n2dhits_hh", ele2dHits, ele_trk->getChi2Ndf(), weight);
595
596 _vtx_histos->Fill2DHisto("pos_track_chi2ndf_v_time_hh", pos_trk->getTrackTime(), pos_trk->getChi2Ndf(), weight);
597 _vtx_histos->Fill2DHisto("pos_track_chi2ndf_v_p_hh", pos_trk->getP(), pos_trk->getChi2Ndf(), weight);
598 _vtx_histos->Fill2DHisto("pos_track_chi2ndf_v_tanlambda_hh", pos_trk->getTanLambda(), pos_trk->getChi2Ndf(), weight);
599 _vtx_histos->Fill2DHisto("pos_track_chi2ndf_v_n2dhits_hh", pos2dHits, pos_trk->getChi2Ndf(), weight);
600
601 _vtx_histos->Fill2DHisto("ele_track_p_v_tanlambda_hh", ele_trk->getTanLambda(), ele_trk->getP(), weight);
602 _vtx_histos->Fill2DHisto("pos_track_p_v_tanlambda_hh", pos_trk->getTanLambda(), pos_trk->getP(), weight);
603
604
605 //1d histos
606 _vtx_histos->Fill1DHisto("ele_track_clus_dt_h", ele_trk->getTrackTime() - corr_eleClusterTime, weight);
607 _vtx_histos->Fill1DHisto("pos_track_clus_dt_h", pos_trk->getTrackTime() - corr_posClusterTime, weight);
608
609 passVtxPresel = true;
610
611 selected_vtxs.push_back(vtx);
612 vtxSelector->clearSelector();
613 }
614
615 // std::cout << "Number of selected vtxs: " << selected_vtxs.size() << std::endl;
616
617 _vtx_histos->Fill1DHisto("n_vertices_h",selected_vtxs.size());
618 if (trks_)
619 _vtx_histos->Fill1DHisto("n_tracks_h",trks_->size());
620
621
622 //not working atm
623 //hps_evt->addVertexCollection("selected_vtxs", selected_vtxs);
624
625 //Make Plots for each region: loop on each region. Check if the region has the cut and apply it
626 //TODO Clean this up => Cuts should be implemented in each region?
627 //TODO Bring the preselection out of this stupid loop
628
629
630 //TODO add yields. => Quite terrible way to loop.
631 for (auto region : _regions ) {
632
633 int nGoodVtx = 0;
634 Vertex* goodVtx = nullptr;
635 std::vector<Vertex*> goodVtxs;
636
637 float truePsum = -1;
638 float trueEsum = -1;
639
640 for ( auto vtx : selected_vtxs) {
641
642 //No cuts.
643 _reg_vtx_selectors[region]->getCutFlowHisto()->Fill(0.,weight);
644
645
646 Particle* ele = nullptr;
647 Particle* pos = nullptr;
648
649 _ah->GetParticlesFromVtx(vtx,ele,pos);
650
651 CalCluster eleClus = ele->getCluster();
652 CalCluster posClus = pos->getCluster();
653 //vtx Z position
654 if (!_reg_vtx_selectors[region]->passCutGt("uncVtxZ_gt",vtx->getZ(),weight))
655 continue;
656
657 double ele_E = ele->getEnergy();
658 double pos_E = pos->getEnergy();
659
660 //Compute analysis variables here.
661
662 Track ele_trk = ele->getTrack();
663 Track pos_trk = pos->getTrack();
664
665 //Beam Position Corrections
666 ele_trk.applyCorrection("z0", beamPosCorrections_.at(1));
667 pos_trk.applyCorrection("z0", beamPosCorrections_.at(1));
668 //Track Time Corrections
669 ele_trk.applyCorrection("track_time",eleTrackTimeBias_);
670 pos_trk.applyCorrection("track_time", posTrackTimeBias_);
671
672 //Get the shared info - TODO change and improve
673
674 Track* ele_trk_gbl = nullptr;
675 Track* pos_trk_gbl = nullptr;
676
677 if (!trkColl_.empty()) {
678 bool foundTracks = _ah->MatchToGBLTracks(ele_trk.getID(),pos_trk.getID(),
679 ele_trk_gbl, pos_trk_gbl, *trks_);
680
681 if (!foundTracks) {
682 if (debug_)
683 std::cout<<"VertexAnaProcessor::ERROR couldn't find ele/pos in the "<<trkColl_ <<"collection"<<std::endl;
684 continue;
685 }
686 }
687 else {
688
689 ele_trk_gbl = (Track*) ele_trk.Clone();
690 pos_trk_gbl = (Track*) pos_trk.Clone();
691 }
692
693 //Add the momenta to the tracks
694 //ele_trk_gbl->setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]);
695 //pos_trk_gbl->setMomentum(pos->getMomentum()[0],pos->getMomentum()[1],pos->getMomentum()[2]);
696 TVector3 recEleP(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]);
697 TLorentzVector p_ele;
698 p_ele.SetPxPyPzE(ele_trk_gbl->getMomentum()[0],ele_trk_gbl->getMomentum()[1],ele_trk_gbl->getMomentum()[2], ele_E);
699 TLorentzVector p_pos;
700 p_pos.SetPxPyPzE(pos_trk_gbl->getMomentum()[0],pos_trk_gbl->getMomentum()[1],pos_trk_gbl->getMomentum()[2], pos_E);
701
702 //Defining these here so they are in scope elsewhere
703 TVector3 trueEleP;
704 TVector3 truePosP;
705
706 if (debug_) {
707 std::cout<<"Check on ele_Track"<<std::endl;
708 std::cout<<"Number of hits:"<<ele_trk_gbl->getTrackerHitCount()<<std::endl;
709 }
710
711 bool foundL1ele = false;
712 bool foundL2ele = false;
713 _ah->InnermostLayerCheck(ele_trk_gbl, foundL1ele, foundL2ele);
714
715
716 if (debug_) {
717 std::cout<<"Check on pos_Track"<<std::endl;
718 std::cout<<"Number of hits:"<<ele_trk_gbl->getTrackerHitCount()<<std::endl;
719 }
720 bool foundL1pos = false;
721 bool foundL2pos = false;
722
723 _ah->InnermostLayerCheck(pos_trk_gbl, foundL1pos, foundL2pos);
724
725 if (debug_) {
726 std::cout<<"Check on pos_Track"<<std::endl;
727 std::cout<<"Innermost:"<<foundL1pos<<" Second Innermost:"<<foundL2pos<<std::endl;
728 }
729
730 //PRESELECTION CUTS
731 if (isData_) {
732 if (!_reg_vtx_selectors[region]->passCutEq("Pair1_eq",(int)evth_->isPair1Trigger(),weight))
733 break;
734 }
735 //Ele Track Time
736 if (!_reg_vtx_selectors[region]->passCutLt("eleTrkTime_lt",fabs(ele_trk_gbl->getTrackTime()),weight))
737 continue;
738
739 //Pos Track Time
740 if (!_reg_vtx_selectors[region]->passCutLt("posTrkTime_lt",fabs(pos_trk_gbl->getTrackTime()),weight))
741 continue;
742
743 //Ele Track-cluster match
744 if (!_reg_vtx_selectors[region]->passCutLt("eleTrkCluMatch_lt",ele->getGoodnessOfPID(),weight))
745 continue;
746
747 //Pos Track-cluster match
748 if (!_reg_vtx_selectors[region]->passCutLt("posTrkCluMatch_lt",pos->getGoodnessOfPID(),weight))
749 continue;
750
751 //Require Positron Cluster exists
752 if (!_reg_vtx_selectors[region]->passCutGt("posClusE_gt",posClus.getEnergy(),weight))
753 continue;
754
755 //Require Positron Cluster does NOT exists
756 if (!_reg_vtx_selectors[region]->passCutLt("posClusE_lt",posClus.getEnergy(),weight))
757 continue;
758
759 double corr_eleClusterTime = ele->getCluster().getTime() - timeOffset_;
760 double corr_posClusterTime = pos->getCluster().getTime() - timeOffset_;
761
762 double botClusTime = 0.0;
763 if(ele->getCluster().getPosition().at(1) < 0.0) botClusTime = ele->getCluster().getTime();
764 else botClusTime = pos->getCluster().getTime();
765
766 //Bottom Cluster Time
767 if (!_reg_vtx_selectors[region]->passCutLt("botCluTime_lt", botClusTime, weight))
768 continue;
769
770 if (!_reg_vtx_selectors[region]->passCutGt("botCluTime_gt", botClusTime, weight))
771 continue;
772
773 //Ele Pos Cluster Time Difference
774 if (!_reg_vtx_selectors[region]->passCutLt("eleposCluTimeDiff_lt",fabs(corr_eleClusterTime - corr_posClusterTime),weight))
775 continue;
776
777 //Ele Track-Cluster Time Difference
778 if (!_reg_vtx_selectors[region]->passCutLt("eleTrkCluTimeDiff_lt",fabs(ele_trk_gbl->getTrackTime() - corr_eleClusterTime),weight))
779 continue;
780
781 //Pos Track-Cluster Time Difference
782 if (!_reg_vtx_selectors[region]->passCutLt("posTrkCluTimeDiff_lt",fabs(pos_trk_gbl->getTrackTime() - corr_posClusterTime),weight))
783 continue;
784
785 TVector3 ele_mom;
786 ele_mom.SetX(ele_trk_gbl->getMomentum()[0]);
787 ele_mom.SetY(ele_trk_gbl->getMomentum()[1]);
788 ele_mom.SetZ(ele_trk_gbl->getMomentum()[2]);
789
790
791 TVector3 pos_mom;
792 pos_mom.SetX(pos_trk_gbl->getMomentum()[0]);
793 pos_mom.SetY(pos_trk_gbl->getMomentum()[1]);
794 pos_mom.SetZ(pos_trk_gbl->getMomentum()[2]);
795
796 //Ele Track Quality - Chi2
797 if (!_reg_vtx_selectors[region]->passCutLt("eleTrkChi2_lt",ele_trk_gbl->getChi2(),weight))
798 continue;
799
800 //Pos Track Quality - Chi2
801 if (!_reg_vtx_selectors[region]->passCutLt("posTrkChi2_lt",pos_trk_gbl->getChi2(),weight))
802 continue;
803
804 //Ele Track Quality - Chi2Ndf
805 if (!_reg_vtx_selectors[region]->passCutLt("eleTrkChi2Ndf_lt",ele_trk_gbl->getChi2Ndf(),weight))
806 continue;
807
808 //Pos Track Quality - Chi2Ndf
809 if (!_reg_vtx_selectors[region]->passCutLt("posTrkChi2Ndf_lt",pos_trk_gbl->getChi2Ndf(),weight))
810 continue;
811
812 //Beam Electron cut
813 if (!_reg_vtx_selectors[region]->passCutLt("eleMom_lt",ele_mom.Mag(),weight))
814 continue;
815
816 //Ele min momentum cut
817 if (!_reg_vtx_selectors[region]->passCutGt("eleMom_gt",ele_mom.Mag(),weight))
818 continue;
819
820 //Pos min momentum cut
821 if (!_reg_vtx_selectors[region]->passCutGt("posMom_gt",pos_mom.Mag(),weight))
822 continue;
823
824 //Ele nHits
825 int ele2dHits = ele_trk_gbl->getTrackerHitCount();
826 if (!ele_trk_gbl->isKalmanTrack())
827 ele2dHits*=2;
828
829 if (!_reg_vtx_selectors[region]->passCutGt("eleN2Dhits_gt",ele2dHits,weight)) {
830 continue;
831 }
832
833 //Pos nHits
834 int pos2dHits = pos_trk_gbl->getTrackerHitCount();
835 if (!pos_trk_gbl->isKalmanTrack())
836 pos2dHits*=2;
837
838 if (!_reg_vtx_selectors[region]->passCutGt("posN2Dhits_gt",pos2dHits,weight)) {
839 continue;
840 }
841
842 //Less than 4 shared hits for ele/pos track
843 if (!_reg_vtx_selectors[region]->passCutLt("eleNshared_lt",ele_trk_gbl->getNShared(),weight)) {
844 continue;
845 }
846
847 if (!_reg_vtx_selectors[region]->passCutLt("posNshared_lt",pos_trk_gbl->getNShared(),weight)) {
848 continue;
849 }
850
851 //Vertex Quality
852 if (!_reg_vtx_selectors[region]->passCutLt("chi2unc_lt",vtx->getChi2(),weight))
853 continue;
854
855 //Max vtx momentum
856 if (!_reg_vtx_selectors[region]->passCutLt("maxVtxMom_lt",(ele_mom+pos_mom).Mag(),weight))
857 continue;
858
859 //Min vtx momentum
860 if (!_reg_vtx_selectors[region]->passCutGt("minVtxMom_gt",(ele_mom+pos_mom).Mag(),weight))
861 continue;
862
863 //END PRESELECTION CUTS
864
865 //L1 requirement
866 if (!_reg_vtx_selectors[region]->passCutEq("L1Requirement_eq",(int)(foundL1ele&&foundL1pos),weight))
867 continue;
868
869 //L2 requirement
870 if (!_reg_vtx_selectors[region]->passCutEq("L2Requirement_eq",(int)(foundL2ele&&foundL2pos),weight))
871 continue;
872
873 //L1 requirement for positron
874 if (!_reg_vtx_selectors[region]->passCutEq("L1PosReq_eq",(int)(foundL1pos),weight))
875 continue;
876
877 //ESum low cut
878 if (!_reg_vtx_selectors[region]->passCutLt("eSum_lt",(ele_E+pos_E),weight))
879 continue;
880
881 //ESum high cut
882 if (!_reg_vtx_selectors[region]->passCutGt("eSum_gt",(ele_E+pos_E),weight))
883 continue;
884
885 //PSum low cut
886 if (!_reg_vtx_selectors[region]->passCutLt("pSum_lt",(p_ele.P()+p_pos.P()),weight))
887 continue;
888
889 //PSum high cut
890 if (!_reg_vtx_selectors[region]->passCutGt("pSum_gt",(p_ele.P()+p_pos.P()),weight))
891 continue;
892
893 //Require Electron Cluster exists
894 if (!_reg_vtx_selectors[region]->passCutGt("eleClusE_gt",eleClus.getEnergy(),weight))
895 continue;
896
897
898 //Require Electron Cluster does NOT exists
899 if (!_reg_vtx_selectors[region]->passCutLt("eleClusE_lt",eleClus.getEnergy(),weight))
900 continue;
901
902 //No shared hits requirement
903 if (!_reg_vtx_selectors[region]->passCutEq("ele_sharedL0_eq",(int)ele_trk_gbl->getSharedLy0(),weight))
904 continue;
905 if (!_reg_vtx_selectors[region]->passCutEq("pos_sharedL0_eq",(int)pos_trk_gbl->getSharedLy0(),weight))
906 continue;
907 if (!_reg_vtx_selectors[region]->passCutEq("ele_sharedL1_eq",(int)ele_trk_gbl->getSharedLy1(),weight))
908 continue;
909 if (!_reg_vtx_selectors[region]->passCutEq("pos_sharedL1_eq",(int)pos_trk_gbl->getSharedLy1(),weight))
910 continue;
911
912 //Min vtx Y pos
913 if (!_reg_vtx_selectors[region]->passCutGt("VtxYPos_gt", vtx->getY(), weight))
914 continue;
915
916 //Max vtx Y pos
917 if (!_reg_vtx_selectors[region]->passCutLt("VtxYPos_lt", vtx->getY(), weight))
918 continue;
919
920 //Tracking Volume for positron
921 if (!_reg_vtx_selectors[region]->passCutGt("volPos_top", p_pos.Py(), weight))
922 continue;
923
924 if (!_reg_vtx_selectors[region]->passCutLt("volPos_bot", p_pos.Py(), weight))
925 continue;
926
927 //If this is MC check if MCParticle matched to the electron track is from rad or recoil
928 if(!isData_)
929 {
930
931 //Fill MC plots after all selections
932 if (!isData_ && mc_reg_on_) _reg_mc_vtx_histos[region]->FillMCParticles(mcParts_, analysis_);
933
934 //Build map of hits and the associated MC part ids for later
935 TRefArray ele_trk_hits = ele_trk_gbl->getSvtHits();
936 TRefArray pos_trk_hits = pos_trk_gbl->getSvtHits();
937 std::map<int, std::vector<int> > trueHitIDs;
938 for(int i = 0; i < hits_->size(); i++)
939 {
940 TrackerHit* hit = hits_->at(i);
941 trueHitIDs[hit->getID()] = hit->getMCPartIDs();
942 }
943 //Count the number of hits per part on the ele track
944 std::map<int, int> nHits4part;
945 for(int i = 0; i < ele_trk_hits.GetEntries(); i++)
946 {
947 TrackerHit* eleHit = (TrackerHit*)ele_trk_hits.At(i);
948 for(int idI = 0; idI < trueHitIDs[eleHit->getID()].size(); idI++ )
949 {
950 int partID = trueHitIDs[eleHit->getID()].at(idI);
951 if ( nHits4part.find(partID) == nHits4part.end() )
952 {
953 // not found
954 nHits4part[partID] = 1;
955 }
956 else
957 {
958 // found
959 nHits4part[partID]++;
960 }
961 }
962 }
963
964 //Determine the MC part with the most hits on the track
965 int maxNHits = 0;
966 int maxID = 0;
967 for (std::map<int,int>::iterator it=nHits4part.begin(); it!=nHits4part.end(); ++it)
968 {
969 if(it->second > maxNHits)
970 {
971 maxNHits = it->second;
972 maxID = it->first;
973 }
974 }
975
976 //Find the correct mc part and grab mother id
977 int isRadEle = -999;
978 int isRecEle = -999;
979
980
981 trueEleP.SetXYZ(-999,-999,-999);
982 truePosP.SetXYZ(-999,-999,-999);
983 if (mcParts_) {
984 float trueEleE = -1;
985 float truePosE = -1;
986 for(int i = 0; i < mcParts_->size(); i++)
987 {
988 int momPDG = mcParts_->at(i)->getMomPDG();
989 if(mcParts_->at(i)->getPDG() == 11 && momPDG == isRadPDG_)
990 {
991 std::vector<double> lP = mcParts_->at(i)->getMomentum();
992 trueEleP.SetXYZ(lP[0],lP[1],lP[2]);
993 trueEleE = mcParts_->at(i)->getEnergy();
994
995 }
996 if(mcParts_->at(i)->getPDG() == -11 && momPDG == isRadPDG_)
997 {
998 std::vector<double> lP = mcParts_->at(i)->getMomentum();
999 truePosP.SetXYZ(lP[0],lP[1],lP[2]);
1000 truePosE = mcParts_->at(i)->getEnergy();
1001
1002 }
1003 if(trueEleP.X() != -999 && truePosP.X() != -999){
1004 truePsum = trueEleP.Mag() + trueEleP.Mag();
1005 trueEsum = trueEleE + truePosE;
1006 }
1007
1008 if(mcParts_->at(i)->getID() != maxID) continue;
1009 //Default isRadPDG = 622
1010 if(momPDG == isRadPDG_) isRadEle = 1;
1011 if(momPDG == 623) isRecEle = 1;
1012 }
1013 }
1014 double momRatio = recEleP.Mag() / trueEleP.Mag();
1015 double momAngle = trueEleP.Angle(recEleP) * TMath::RadToDeg();
1016 if (!_reg_vtx_selectors[region]->passCutLt("momRatio_lt", momRatio, weight)) continue;
1017 if (!_reg_vtx_selectors[region]->passCutGt("momRatio_gt", momRatio, weight)) continue;
1018 if (!_reg_vtx_selectors[region]->passCutLt("momAngle_lt", momAngle, weight)) continue;
1019
1020 if (!_reg_vtx_selectors[region]->passCutEq("isRadEle_eq", isRadEle, weight)) continue;
1021 if (!_reg_vtx_selectors[region]->passCutEq("isNotRadEle_eq", isRadEle, weight)) continue;
1022 if (!_reg_vtx_selectors[region]->passCutEq("isRecEle_eq", isRecEle, weight)) continue;
1023 }
1024
1025 goodVtx = vtx;
1026 nGoodVtx++;
1027 goodVtxs.push_back(vtx);
1028 } // selected vertices
1029
1030 //N selected vertices - this is quite a silly cut to make at the end. But okay. that's how we decided atm.
1031 if (!_reg_vtx_selectors[region]->passCutEq("nVtxs_eq", nGoodVtx, weight))
1032 continue;
1033 //Move to after N vertices cut (was filled before)
1034 _reg_vtx_histos[region]->Fill1DHisto("n_vertices_h", nGoodVtx, weight);
1035
1036 //Loop over all selected vertices in the region
1037 for(std::vector<Vertex*>::iterator it = goodVtxs.begin(); it != goodVtxs.end(); it++){
1038
1039 Vertex* vtx = *it;
1040
1041 Particle* ele = nullptr;
1042 Particle* pos = nullptr;
1043
1044 if (!vtx || !_ah->GetParticlesFromVtx(vtx,ele,pos))
1045 continue;
1046
1047 CalCluster eleClus = ele->getCluster();
1048 CalCluster posClus = pos->getCluster();
1049
1050 double corr_eleClusterTime = ele->getCluster().getTime() - timeOffset_;
1051 double corr_posClusterTime = pos->getCluster().getTime() - timeOffset_;
1052
1053 double ele_E = ele->getEnergy();
1054 double pos_E = pos->getEnergy();
1055
1056 //Compute analysis variables here.
1057 Track ele_trk = ele->getTrack();
1058 Track pos_trk = pos->getTrack();
1059 //Get the shared info - TODO change and improve
1060
1061 Track* ele_trk_gbl = nullptr;
1062 Track* pos_trk_gbl = nullptr;
1063
1064 if (!trkColl_.empty()) {
1065 bool foundTracks = _ah->MatchToGBLTracks(ele_trk.getID(),pos_trk.getID(),
1066 ele_trk_gbl, pos_trk_gbl, *trks_);
1067
1068 if (!foundTracks) {
1069 if (debug_)
1070 std::cout<<"VertexAnaProcessor::ERROR couldn't find ele/pos in the "<<trkColl_ <<"collection"<<std::endl;
1071 continue;
1072 }
1073 }
1074 else {
1075
1076 ele_trk_gbl = (Track*) ele_trk.Clone();
1077 pos_trk_gbl = (Track*) pos_trk.Clone();
1078 }
1079
1080 //Vertex Covariance
1081 std::vector<float> vtx_cov = vtx->getCovariance();
1082 float cxx = vtx_cov.at(0);
1083 float cyx = vtx_cov.at(1);
1084 float cyy = vtx_cov.at(2);
1085 float czx = vtx_cov.at(3);
1086 float czy = vtx_cov.at(4);
1087 float czz = vtx_cov.at(5);
1088
1089 //MC Truth hits in first 4 sensors
1090 int L1L2hitCode = 0; //hit code '1111' means truth ax+ster hits in L1_ele, L1_pos, L2_ele, L2_pos
1091 int L1hitCode = 0; //hit code '1111' means truth in L1_ele_ax, L1_ele_ster, L1_pos_ax, L1_pos_ster
1092 int L2hitCode = 0; // hit code '1111' means truth in L2_ele_ax, L2_ele_ster, L2_pos_ax, L2_pos_ster
1093 if(!isData_){
1094 //Get hit codes. Only sure this works for 2016 KF as is.
1095 utils::get2016KFMCTruthHitCodes(ele_trk_gbl, pos_trk_gbl, L1L2hitCode, L1hitCode, L2hitCode);
1096 //L1L2 truth hit selection
1097 if (!_reg_vtx_selectors[region]->passCutLt("hitCode_lt",((double)L1L2hitCode)-0.5, weight)) continue;
1098 if (!_reg_vtx_selectors[region]->passCutGt("hitCode_gt",((double)L1L2hitCode)+0.5, weight)) continue;
1099 //Fil hitcodes
1100 _reg_vtx_histos[region]->Fill1DHisto("hitCode_h", L1L2hitCode,weight);
1101 _reg_vtx_histos[region]->Fill1DHisto("L1hitCode_h", L1hitCode,weight);
1102 _reg_vtx_histos[region]->Fill1DHisto("L2hitCode_h", L2hitCode,weight);
1103 }
1104
1105 //track isolations
1106 //Only calculate isolations if both track L1 and L2 hits exist
1107 bool hasL1ele = false;
1108 bool hasL2ele = false;
1109 _ah->InnermostLayerCheck(ele_trk_gbl, hasL1ele, hasL2ele);
1110
1111 bool hasL1pos = false;
1112 bool hasL2pos = false;
1113 _ah->InnermostLayerCheck(pos_trk_gbl, hasL1pos, hasL2pos);
1114
1115 double ele_trk_iso_L1 = 99999.9;
1116 double pos_trk_iso_L1 = 99999.9;
1117 if(hasL1ele && hasL2ele && hasL1pos && hasL2pos){
1118 if (ele_trk_gbl->isKalmanTrack()){
1119 ele_trk_iso_L1 = utils::getKalmanTrackL1Isolations(ele_trk_gbl, hits_);
1120 pos_trk_iso_L1 = utils::getKalmanTrackL1Isolations(pos_trk_gbl, hits_);
1121 }
1122 }
1123
1124 TVector3 ele_mom;
1125 //ele_mom.SetX(ele->getMomentum()[0]);
1126 //ele_mom.SetY(ele->getMomentum()[1]);
1127 //ele_mom.SetZ(ele->getMomentum()[2]);
1128 ele_mom.SetX(ele_trk_gbl->getMomentum()[0]);
1129 ele_mom.SetY(ele_trk_gbl->getMomentum()[1]);
1130 ele_mom.SetZ(ele_trk_gbl->getMomentum()[2]);
1131
1132
1133 TVector3 pos_mom;
1134 //pos_mom.SetX(pos->getMomentum()[0]);
1135 //pos_mom.SetY(pos->getMomentum()[1]);
1136 //pos_mom.SetZ(pos->getMomentum()[2]);
1137 pos_mom.SetX(pos_trk_gbl->getMomentum()[0]);
1138 pos_mom.SetY(pos_trk_gbl->getMomentum()[1]);
1139 pos_mom.SetZ(pos_trk_gbl->getMomentum()[2]);
1140
1141 double ele_pos_dt = corr_eleClusterTime - corr_posClusterTime;
1142 double psum = ele_mom.Mag()+pos_mom.Mag();
1143
1144 //Ele nHits
1145 int ele2dHits = ele_trk_gbl->getTrackerHitCount();
1146 if (!ele_trk_gbl->isKalmanTrack())
1147 ele2dHits*=2;
1148
1149 //pos nHits
1150 int pos2dHits = pos_trk_gbl->getTrackerHitCount();
1151 if (!pos_trk_gbl->isKalmanTrack())
1152 pos2dHits*=2;
1153
1154 if(ts_ != nullptr)
1155 {
1156 _reg_vtx_histos[region]->Fill2DHisto("trig_count_hh",
1157 ((int)ts_->prescaled.Single_3_Top)+((int)ts_->prescaled.Single_3_Bot),
1158 ((int)ts_->prescaled.Single_2_Top)+((int)ts_->prescaled.Single_2_Bot));
1159 }
1160 _reg_vtx_histos[region]->Fill1DHisto("n_vtx_h", vtxs_->size());
1161 int NposTrks = 0;
1162 int NeleTrks = 0;
1163 for (int iT = 0; iT < trks_->size(); iT++)
1164 {
1165 if (trks_->at(iT)->getCharge() > 0) NposTrks++;
1166 else NeleTrks++;
1167 }
1168 _reg_vtx_histos[region]->Fill2DHisto("n_tracks_hh", NeleTrks, NposTrks);
1169
1170 //Add the momenta to the tracks
1171 //ele_trk_gbl->setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]);
1172 //pos_trk_gbl->setMomentum(pos->getMomentum()[0],pos->getMomentum()[1],pos->getMomentum()[2]);
1173 TVector3 recEleP(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]);
1174 TLorentzVector p_ele;
1175 p_ele.SetPxPyPzE(ele_trk_gbl->getMomentum()[0],ele_trk_gbl->getMomentum()[1],ele_trk_gbl->getMomentum()[2], ele_E);
1176 TLorentzVector p_pos;
1177 p_pos.SetPxPyPzE(pos_trk_gbl->getMomentum()[0],pos_trk_gbl->getMomentum()[1],pos_trk_gbl->getMomentum()[2], pos_E);
1178
1179
1180 _reg_vtx_histos[region]->Fill2DHistograms(vtx,weight);
1181 _reg_vtx_histos[region]->Fill1DVertex(vtx,
1182 ele,
1183 pos,
1184 ele_trk_gbl,
1185 pos_trk_gbl,
1186 weight);
1187
1188 _reg_vtx_histos[region]->Fill1DHisto("ele_pos_clusTimeDiff_h", (corr_eleClusterTime - corr_posClusterTime), weight);
1189 _reg_vtx_histos[region]->Fill1DHisto("ele_track_n2dhits_h", ele2dHits, weight);
1190 _reg_vtx_histos[region]->Fill1DHisto("pos_track_n2dhits_h", pos2dHits, weight);
1191 _reg_vtx_histos[region]->Fill1DHisto("vtx_Psum_h", p_ele.P()+p_pos.P(), weight);
1192 _reg_vtx_histos[region]->Fill1DHisto("vtx_Esum_h", eleClus.getEnergy()+posClus.getEnergy(), weight);
1193 _reg_vtx_histos[region]->Fill2DHisto("ele_vtxZ_iso_hh", TMath::Min(ele_trk_gbl->getIsolation(0), ele_trk_gbl->getIsolation(1)), vtx->getZ(), weight);
1194 _reg_vtx_histos[region]->Fill2DHisto("pos_vtxZ_iso_hh", TMath::Min(pos_trk_gbl->getIsolation(0), pos_trk_gbl->getIsolation(1)), vtx->getZ(), weight);
1195 _reg_vtx_histos[region]->Fill2DTrack(ele_trk_gbl,weight,"ele_");
1196 _reg_vtx_histos[region]->Fill2DTrack(pos_trk_gbl,weight,"pos_");
1197 _reg_vtx_histos[region]->Fill1DHisto("mcMass622_h",apMass);
1198 _reg_vtx_histos[region]->Fill1DHisto("mcZ622_h",apZ);
1199 _reg_vtx_histos[region]->Fill1DHisto("mcMass625_h",vdMass);
1200 _reg_vtx_histos[region]->Fill1DHisto("mcZ625_h",vdZ);
1201
1202 if (trks_) _reg_vtx_histos[region]->Fill1DHisto("n_tracks_h",trks_->size(),weight);
1203 _reg_vtx_histos[region]->Fill1DHisto("vtx_track_L1_isolation_h", ele_trk_iso_L1);
1204 _reg_vtx_histos[region]->Fill1DHisto("vtx_track_L1_isolation_h", pos_trk_iso_L1);
1205
1206 //Just for the selected vertex
1207 if(!isData_)
1208 {
1209 _reg_vtx_histos[region]->Fill2DHisto("vtx_Esum_vs_true_Esum_hh",eleClus.getEnergy()+posClus.getEnergy(), trueEsum, weight);
1210 _reg_vtx_histos[region]->Fill2DHisto("vtx_Psum_vs_true_Psum_hh",p_ele.P()+p_pos.P(), truePsum, weight);
1211 _reg_vtx_histos[region]->Fill1DHisto("true_vtx_psum_h",truePsum,weight);
1212 }
1213
1214 double reconz = vtx->getZ();
1215 double ele_trk_z0 = ele_trk_gbl->getZ0();
1216 double ele_trk_z0err = ele_trk_gbl->getZ0Err();
1217 double pos_trk_z0 = pos_trk_gbl->getZ0();
1218 double pos_trk_z0err = pos_trk_gbl->getZ0Err();
1219
1220 //Project vertex to target
1221 double vtx_proj_x = -999.9;
1222 double vtx_proj_y = -999.9;
1223 double vtx_proj_x_sig = -999.9;
1224 double vtx_proj_y_sig = -999.9;
1225 double vtx_proj_sig = -999.9;
1226 if(!v0ProjectionFitsCfg_.empty())
1227 vtx_proj_sig = utils::v0_projection_to_target_significance(v0proj_fits_, evth_->getRunNumber(),
1228 vtx_proj_x, vtx_proj_y, vtx_proj_x_sig, vtx_proj_y_sig, vtx->getX(), vtx->getY(),
1229 reconz, vtx->getP().X(), vtx->getP().Y(), vtx->getP().Z());
1230
1231 _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_x_v_unc_vtx_y_hh", vtx->getX(), vtx->getY());
1232 _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_proj_x_v_unc_vtx_proj_y_hh", vtx_proj_x, vtx_proj_y);
1233 _reg_vtx_histos[region]->Fill2DHisto("unc_vtx_proj_x_y_significance_hh", vtx_proj_x_sig, vtx_proj_y_sig);
1234 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_vtx_proj_significance_hh", vtx_proj_sig, reconz);
1235 _reg_vtx_histos[region]->Fill2DHisto("vtx_track_reconz_v_Z0err_hh", ele_trk_z0err, reconz);
1236 _reg_vtx_histos[region]->Fill2DHisto("vtx_track_reconz_v_Z0err_hh", pos_trk_z0err, reconz);
1237 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_z0_hh", ele_trk_z0, reconz);
1238 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_z0_hh", pos_trk_z0, reconz);
1239 _reg_vtx_histos[region]->Fill2DHisto("vtx_track_recon_z_v_ABSdz0tanlambda_hh", std::abs((ele_trk_z0/ele_trk_gbl->getTanLambda()) - (pos_trk_z0/pos_trk_gbl->getTanLambda())), reconz);
1240 _reg_vtx_histos[region]->Fill2DHisto("vtx_track_recon_z_v_dz0tanlambda_hh", ((ele_trk_z0/ele_trk_gbl->getTanLambda()) - (pos_trk_z0/pos_trk_gbl->getTanLambda())), reconz);
1241 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_cxx_hh", cxx, reconz);
1242 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_cyy_hh", cyy, reconz);
1243 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_czz_hh", czz, reconz);
1244 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_cyx_hh", cyx, reconz);
1245 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_czx_hh", czx, reconz);
1246 _reg_vtx_histos[region]->Fill2DHisto("recon_z_v_czy_hh", czy, reconz);
1247 _reg_vtx_histos[region]->Fill1DHisto("cxx_h", cxx);
1248 _reg_vtx_histos[region]->Fill1DHisto("cyy_h", cyy);
1249 _reg_vtx_histos[region]->Fill1DHisto("czz_h", czz);
1250 _reg_vtx_histos[region]->Fill1DHisto("cyx_h", cyx);
1251 _reg_vtx_histos[region]->Fill1DHisto("czx_h", czx);
1252 _reg_vtx_histos[region]->Fill1DHisto("czy_h", czy);
1253 _reg_vtx_histos[region]->Fill2DHisto("vtx_track_recon_z_v_z0tanlambda_hh", ele_trk_z0/ele_trk_gbl->getTanLambda(), reconz);
1254 _reg_vtx_histos[region]->Fill2DHisto("vtx_track_recon_z_v_z0tanlambda_hh", pos_trk_z0/pos_trk_gbl->getTanLambda(), reconz);
1255 _reg_vtx_histos[region]->Fill2DHisto("ele_clusT_v_ele_trackT_hh", ele_trk_gbl->getTrackTime(), corr_eleClusterTime, weight);
1256 _reg_vtx_histos[region]->Fill2DHisto("pos_clusT_v_pos_trackT_hh", pos_trk_gbl->getTrackTime(), corr_posClusterTime, weight);
1257 _reg_vtx_histos[region]->Fill2DHisto("ele_track_time_v_P_hh", ele_trk_gbl->getP(), ele_trk_gbl->getTrackTime(), weight);
1258 _reg_vtx_histos[region]->Fill2DHisto("pos_track_time_v_P_hh", pos_trk_gbl->getP(), pos_trk_gbl->getTrackTime(), weight);
1259 _reg_vtx_histos[region]->Fill2DHisto("ele_pos_clusTimeDiff_v_pSum_hh",ele_mom.Mag()+pos_mom.Mag(), ele_pos_dt, weight);
1260 _reg_vtx_histos[region]->Fill2DHisto("ele_cluster_energy_v_track_p_hh",ele_trk_gbl->getP(), eleClus.getEnergy(), weight);
1261 _reg_vtx_histos[region]->Fill2DHisto("pos_cluster_energy_v_track_p_hh",pos_trk_gbl->getP(), posClus.getEnergy(), weight);
1262 _reg_vtx_histos[region]->Fill2DHisto("ele_track_cluster_dt_v_EoverP_hh",eleClus.getEnergy()/ele_trk_gbl->getP(), ele_trk_gbl->getTrackTime() - corr_eleClusterTime, weight);
1263 _reg_vtx_histos[region]->Fill2DHisto("pos_track_cluster_dt_v_EoverP_hh",posClus.getEnergy()/pos_trk_gbl->getP(), pos_trk_gbl->getTrackTime() - corr_posClusterTime, weight);
1264 _reg_vtx_histos[region]->Fill2DHisto("ele_track_clus_dt_v_p_hh",ele_trk_gbl->getP(), ele_trk_gbl->getTrackTime() - corr_eleClusterTime, weight);
1265 _reg_vtx_histos[region]->Fill2DHisto("pos_track_clus_dt_v_p_hh",pos_trk_gbl->getP(), pos_trk_gbl->getTrackTime() - corr_posClusterTime, weight);
1266 _reg_vtx_histos[region]->Fill2DHisto("ele_z0_vs_pos_z0_hh",ele_trk_gbl->getZ0(), pos_trk_gbl->getZ0(), weight);
1267 //chi2 2d plots
1268 _reg_vtx_histos[region]->Fill2DHisto("ele_track_chi2ndf_v_time_hh", ele_trk_gbl->getTrackTime(), ele_trk_gbl->getChi2Ndf(), weight);
1269 _reg_vtx_histos[region]->Fill2DHisto("ele_track_chi2ndf_v_p_hh", ele_trk_gbl->getP(), ele_trk_gbl->getChi2Ndf(), weight);
1270 _reg_vtx_histos[region]->Fill2DHisto("ele_track_chi2ndf_v_tanlambda_hh", ele_trk_gbl->getTanLambda(), ele_trk_gbl->getChi2Ndf(), weight);
1271 _reg_vtx_histos[region]->Fill2DHisto("ele_track_chi2ndf_v_n2dhits_hh", ele2dHits, ele_trk_gbl->getChi2Ndf(), weight);
1272
1273 _reg_vtx_histos[region]->Fill2DHisto("pos_track_chi2ndf_v_time_hh", pos_trk_gbl->getTrackTime(), pos_trk_gbl->getChi2Ndf(), weight);
1274 _reg_vtx_histos[region]->Fill2DHisto("pos_track_chi2ndf_v_p_hh", pos_trk_gbl->getP(), pos_trk_gbl->getChi2Ndf(), weight);
1275 _reg_vtx_histos[region]->Fill2DHisto("pos_track_chi2ndf_v_tanlambda_hh", pos_trk_gbl->getTanLambda(), pos_trk_gbl->getChi2Ndf(), weight);
1276 _reg_vtx_histos[region]->Fill2DHisto("pos_track_chi2ndf_v_n2dhits_hh", pos2dHits, pos_trk_gbl->getChi2Ndf(), weight);
1277
1278 _reg_vtx_histos[region]->Fill2DHisto("ele_track_p_v_tanlambda_hh", ele_trk_gbl->getTanLambda(), ele_trk_gbl->getP(), weight);
1279 _reg_vtx_histos[region]->Fill2DHisto("pos_track_p_v_tanlambda_hh", pos_trk_gbl->getTanLambda(), pos_trk_gbl->getP(), weight);
1280
1281
1282 //1d histos
1283 _reg_vtx_histos[region]->Fill1DHisto("ele_track_clus_dt_h", ele_trk_gbl->getTrackTime() - corr_eleClusterTime, weight);
1284 _reg_vtx_histos[region]->Fill1DHisto("pos_track_clus_dt_h", pos_trk_gbl->getTrackTime() - corr_posClusterTime, weight);
1285
1286 //TODO put this in the Vertex!
1287 TVector3 vtxPosSvt;
1288 vtxPosSvt.SetX(vtx->getX());
1289 vtxPosSvt.SetY(vtx->getY());
1290 vtxPosSvt.SetZ(vtx->getZ());
1291 vtxPosSvt.RotateY(-0.0305);
1292
1293 //Just for the selected vertex
1294 if (makeFlatTuple_){
1295 if(!isData_){
1296 _reg_tuples[region]->setVariableValue("ap_true_vtx_z", apZ);
1297 _reg_tuples[region]->setVariableValue("ap_true_vtx_mass", apMass);
1298 _reg_tuples[region]->setVariableValue("ap_true_vtx_energy", apEnergy);
1299 _reg_tuples[region]->setVariableValue("vd_true_vtx_z", vdZ);
1300 _reg_tuples[region]->setVariableValue("vd_true_vtx_mass", vdMass);
1301 _reg_tuples[region]->setVariableValue("vd_true_vtx_energy", vdEnergy);
1302 _reg_tuples[region]->setVariableValue("hitCode", float(L1L2hitCode));
1303 _reg_tuples[region]->setVariableValue("L1hitCode", float(L1hitCode));
1304 _reg_tuples[region]->setVariableValue("L2hitCode", float(L2hitCode));
1305 }
1306
1307 _reg_tuples[region]->setVariableValue("unc_vtx_mass", vtx->getInvMass());
1308 _reg_tuples[region]->setVariableValue("unc_vtx_z" , vtxPosSvt.Z());
1309 _reg_tuples[region]->setVariableValue("unc_vtx_chi2", vtx->getChi2());
1310 _reg_tuples[region]->setVariableValue("unc_vtx_psum", p_ele.P()+p_pos.P());
1311 _reg_tuples[region]->setVariableValue("unc_vtx_px", vtx->getP().X());
1312 _reg_tuples[region]->setVariableValue("unc_vtx_py", vtx->getP().Y());
1313 _reg_tuples[region]->setVariableValue("unc_vtx_pz", vtx->getP().Z());
1314 _reg_tuples[region]->setVariableValue("unc_vtx_x", vtx->getX());
1315 _reg_tuples[region]->setVariableValue("unc_vtx_y", vtx->getY());
1316 _reg_tuples[region]->setVariableValue("unc_vtx_ele_pos_clust_dt", corr_eleClusterTime - corr_posClusterTime);
1317 _reg_tuples[region]->setVariableValue("unc_vtx_cxx", cxx);
1318 _reg_tuples[region]->setVariableValue("unc_vtx_cyy", cyy);
1319 _reg_tuples[region]->setVariableValue("unc_vtx_czz", czz);
1320 _reg_tuples[region]->setVariableValue("unc_vtx_cyx", cyx);
1321 _reg_tuples[region]->setVariableValue("unc_vtx_czy", czy);
1322 _reg_tuples[region]->setVariableValue("unc_vtx_czx", czx);
1323 _reg_tuples[region]->setVariableValue("unc_vtx_proj_x", vtx_proj_x);
1324 _reg_tuples[region]->setVariableValue("unc_vtx_proj_y", vtx_proj_y);
1325 _reg_tuples[region]->setVariableValue("unc_vtx_proj_x_sig", vtx_proj_x_sig);
1326 _reg_tuples[region]->setVariableValue("unc_vtx_proj_y_sig", vtx_proj_y_sig);
1327 _reg_tuples[region]->setVariableValue("unc_vtx_proj_sig", vtx_proj_sig);
1328
1329 //track vars
1330 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_p", ele_trk_gbl->getP());
1331 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_t", ele_trk_gbl->getTrackTime());
1332 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_d0", ele_trk_gbl->getD0());
1333 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_phi0", ele_trk_gbl->getPhi());
1334 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_omega", ele_trk_gbl->getOmega());
1335 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_tanLambda", ele_trk_gbl->getTanLambda());
1336 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_z0", ele_trk_gbl->getZ0());
1337 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_chi2ndf", ele_trk_gbl->getChi2Ndf());
1338 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_clust_dt", ele_trk_gbl->getTrackTime() - corr_eleClusterTime);
1339 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_z0Err",ele_trk_gbl->getZ0Err());
1340 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_d0Err", ele_trk_gbl->getD0Err());
1341 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_tanLambdaErr", ele_trk_gbl->getTanLambdaErr());
1342 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_PhiErr", ele_trk_gbl->getPhiErr());
1343 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_OmegaErr", ele_trk_gbl->getOmegaErr());
1344 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_L1_isolation",ele_trk_iso_L1);
1345 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_nhits",ele2dHits);
1346
1347 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_p", pos_trk_gbl->getP());
1348 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_t", pos_trk_gbl->getTrackTime());
1349 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_d0", pos_trk_gbl->getD0());
1350 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_phi0", pos_trk_gbl->getPhi());
1351 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_omega", pos_trk_gbl->getOmega());
1352 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_tanLambda", pos_trk_gbl->getTanLambda());
1353 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_z0", pos_trk_gbl->getZ0());
1354 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_chi2ndf", pos_trk_gbl->getChi2Ndf());
1355 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_clust_dt", pos_trk_gbl->getTrackTime() - corr_posClusterTime);
1356 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_z0Err",pos_trk_gbl->getZ0Err());
1357 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_d0Err", pos_trk_gbl->getD0Err());
1358 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_tanLambdaErr", pos_trk_gbl->getTanLambdaErr());
1359 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_PhiErr", pos_trk_gbl->getPhiErr());
1360 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_OmegaErr", pos_trk_gbl->getOmegaErr());
1361 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_L1_isolation",pos_trk_iso_L1);
1362 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_nhits",pos2dHits);
1363
1364 //clust vars
1365 _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_E", eleClus.getEnergy());
1366 _reg_tuples[region]->setVariableValue("unc_vtx_ele_clust_corr_t",corr_eleClusterTime);
1367
1368 _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_E", posClus.getEnergy());
1369 _reg_tuples[region]->setVariableValue("unc_vtx_pos_clust_corr_t",corr_posClusterTime);
1370 _reg_tuples[region]->setVariableValue("run_number", evth_->getRunNumber());
1371
1372 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_x", ele_trk_gbl->getPosition().at(0));
1373 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_y", ele_trk_gbl->getPosition().at(1));
1374 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_z", ele_trk_gbl->getPosition().at(2));
1375 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_x", pos_trk_gbl->getPosition().at(0));
1376 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_y", pos_trk_gbl->getPosition().at(1));
1377 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_z", pos_trk_gbl->getPosition().at(2));
1378 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_px", ele_trk_gbl->getMomentum().at(0));
1379 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_py", ele_trk_gbl->getMomentum().at(1));
1380 _reg_tuples[region]->setVariableValue("unc_vtx_ele_track_pz", ele_trk_gbl->getMomentum().at(2));
1381 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_px", pos_trk_gbl->getMomentum().at(0));
1382 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_py", pos_trk_gbl->getMomentum().at(1));
1383 _reg_tuples[region]->setVariableValue("unc_vtx_pos_track_pz", pos_trk_gbl->getMomentum().at(2));
1384
1385 _reg_tuples[region]->fill();
1386 }
1387 }
1388
1389 }// regions
1390
1391 return true;
1392}
1393
1394void VertexAnaProcessor::finalize() {
1395
1396 //TODO clean this up a little.
1397 outF_->cd();
1398 _vtx_histos->saveHistos(outF_,_vtx_histos->getName());
1399 outF_->cd(_vtx_histos->getName().c_str());
1400 vtxSelector->getCutFlowHisto()->Write();
1401
1402 outF_->cd();
1403 if(!isData_ && mc_reg_on_)
1404 _mc_vtx_histos->saveHistos(outF_, _mc_vtx_histos->getName());
1405 //delete histos;
1406 //histos = nullptr;
1407
1408
1409 for (reg_it it = _reg_vtx_histos.begin(); it!=_reg_vtx_histos.end(); ++it) {
1410 std::string dirName = anaName_+"_"+it->first;
1411 (it->second)->saveHistos(outF_,dirName);
1412 outF_->cd(dirName.c_str());
1413 _reg_vtx_selectors[it->first]->getCutFlowHisto()->Write();
1414 //Save tuples
1415 if (makeFlatTuple_)
1416 _reg_tuples[it->first]->writeTree();
1417
1418 }
1419
1420 if(!isData_ && mc_reg_on_){
1421 for (reg_mc_it it = _reg_mc_vtx_histos.begin(); it!=_reg_mc_vtx_histos.end(); ++it) {
1422 std::string dirName = anaName_+"_mc_"+it->first;
1423 (it->second)->saveHistos(outF_,dirName);
1424 outF_->cd(dirName.c_str());
1425 }
1426 }
1427
1428 outF_->Close();
1429
1430}
1431
1432DECLARE_PROCESSOR(VertexAnaProcessor);
DECLARE_PROCESSOR
#define DECLARE_PROCESSOR(CLASS)
Macro which allows the framework to construct a producer given its name during configuration.
Definition Processor.h:139
AnaHelpers::getFileName
static std::string getFileName(std::string filePath, bool withExtension)
Definition AnaHelpers.cxx:5
CalCluster::getEnergy
double getEnergy() const
Definition CalCluster.h:73
CalCluster::getPosition
std::vector< double > getPosition() const
Definition CalCluster.h:63
CalCluster::getTime
double getTime() const
Definition CalCluster.h:83
EventHeader::isPair1Trigger
bool isPair1Trigger() const
Definition EventHeader.h:101
EventHeader::getEventNumber
int getEventNumber() const
Definition EventHeader.h:50
EventHeader::getRunNumber
int getRunNumber() const
Definition EventHeader.h:69
IEvent
Definition IEvent.h:7
ParameterSet
description
Definition ParameterSet.h:22
Particle::getTrack
Track getTrack() const
Definition Particle.h:48
Particle::getMomentum
std::vector< double > getMomentum() const
Definition Particle.cxx:29
Particle::getCluster
CalCluster getCluster() const
Definition Particle.h:62
Particle::getEnergy
double getEnergy() const
Definition Particle.h:143
Particle::getGoodnessOfPID
double getGoodnessOfPID() const
Definition Particle.h:134
Processor
Base class for all event processing components.
Definition Processor.h:34
Processor::outF_
TFile * outF_
Definition Processor.h:125
Processor::process
virtual bool process()
Process the histograms and generate analysis output.
Definition Processor.h:95
TSData::prescaled
tsBits prescaled
Definition TSData.h:66
Track
Definition Track.h:32
Track::isKalmanTrack
bool isKalmanTrack() const
Definition Track.h:231
Track::getTanLambdaErr
double getTanLambdaErr() const
Definition Track.h:108
Track::getPosition
std::vector< double > getPosition()
Definition Track.cxx:61
Track::getPhiErr
double getPhiErr() const
Definition Track.h:106
Track::getPhi
double getPhi() const
Definition Track.h:88
Track::getChi2Ndf
double getChi2Ndf() const
Definition Track.h:126
Track::getChi2
double getChi2() const
Definition Track.h:123
Track::getTrackerHitCount
int getTrackerHitCount() const
Definition Track.h:345
Track::getOmega
double getOmega() const
Definition Track.h:90
Track::getTrackTime
double getTrackTime() const
Definition Track.h:160
Track::getD0
double getD0() const
Definition Track.h:87
Track::getD0Err
double getD0Err() const
Definition Track.h:105
Track::getSharedLy0
bool getSharedLy0() const
Definition Track.h:355
Track::getID
int getID() const
Definition Track.h:262
Track::getNShared
int getNShared() const
Definition Track.h:350
Track::getTanLambda
double getTanLambda() const
Definition Track.h:91
Track::getSvtHits
TRefArray getSvtHits() const
Definition Track.h:64
Track::getIsolation
double getIsolation(const int layer) const
Definition Track.h:149
Track::getMomentum
std::vector< double > getMomentum()
Definition Track.h:285
Track::getSharedLy1
bool getSharedLy1() const
Definition Track.h:356
Track::getZ0
double getZ0() const
Definition Track.h:92
Track::getP
double getP() const
Definition Track.h:291
Track::getZ0Err
double getZ0Err() const
Definition Track.h:109
Track::applyCorrection
void applyCorrection(std::string var, double correction)
Definition Track.cxx:91
Track::getOmegaErr
double getOmegaErr() const
Definition Track.h:107
TrackerHit::getMCPartIDs
std::vector< int > getMCPartIDs() const
Definition TrackerHit.h:127
TrackerHit::getID
int getID() const
Definition TrackerHit.h:124
VertexAnaProcessor
Insert description here. more details.
Definition VertexAnaProcessor.h:45
VertexAnaProcessor::_mc_vtx_histos
std::shared_ptr< MCAnaHistos > _mc_vtx_histos
description
Definition VertexAnaProcessor.h:121
VertexAnaProcessor::v0ProjectionFitsCfg_
std::string v0ProjectionFitsCfg_
json file w run dependent v0 projection fits
Definition VertexAnaProcessor.h:147
VertexAnaProcessor::bpc_configs_
json bpc_configs_
json object
Definition VertexAnaProcessor.h:145
VertexAnaProcessor::histoCfg_
std::string histoCfg_
description
Definition VertexAnaProcessor.h:134
VertexAnaProcessor::timeOffset_
double timeOffset_
description
Definition VertexAnaProcessor.h:136
VertexAnaProcessor::selectionCfg_
std::string selectionCfg_
description
Definition VertexAnaProcessor.h:91
VertexAnaProcessor::anaName_
std::string anaName_
description
Definition VertexAnaProcessor.h:109
VertexAnaProcessor::makeFlatTuple_
int makeFlatTuple_
make true in config to save flat tuple
Definition VertexAnaProcessor.h:117
VertexAnaProcessor::configure
virtual void configure(const ParameterSet &parameters)
description
Definition VertexAnaProcessor.cxx:20
VertexAnaProcessor::finalize
virtual void finalize()
description
Definition VertexAnaProcessor.cxx:1394
VertexAnaProcessor::ts_
TSData * ts_
description
Definition VertexAnaProcessor.h:102
VertexAnaProcessor::regionSelections_
std::vector< std::string > regionSelections_
description
Definition VertexAnaProcessor.h:89
VertexAnaProcessor::initialize
virtual void initialize(TTree *tree)
description
Definition VertexAnaProcessor.cxx:62
VertexAnaProcessor::_reg_vtx_histos
std::map< std::string, std::shared_ptr< TrackHistos > > _reg_vtx_histos
description
Definition VertexAnaProcessor.h:125
VertexAnaProcessor::reg_mc_it
std::map< std::string, std::shared_ptr< MCAnaHistos > >::iterator reg_mc_it
description
Definition VertexAnaProcessor.h:132
VertexAnaProcessor::analysis_
std::string analysis_
description
Definition VertexAnaProcessor.h:139
VertexAnaProcessor::VertexAnaProcessor
VertexAnaProcessor(const std::string &name, Process &process)
Constructor.
Definition VertexAnaProcessor.cxx:12
VertexAnaProcessor::tree_
TTree * tree_
description
Definition VertexAnaProcessor.h:118
VertexAnaProcessor::bvtxs_
TBranch * bvtxs_
description
Definition VertexAnaProcessor.h:94
VertexAnaProcessor::mcColl_
std::string mcColl_
description
Definition VertexAnaProcessor.h:115
VertexAnaProcessor::becal_
TBranch * becal_
description
Definition VertexAnaProcessor.h:99
VertexAnaProcessor::beamE_
double beamE_
In GeV. Default is 2016 value;.
Definition VertexAnaProcessor.h:137
VertexAnaProcessor::trks_
std::vector< Track * > * trks_
description
Definition VertexAnaProcessor.h:105
VertexAnaProcessor::_reg_vtx_selectors
std::map< std::string, std::shared_ptr< BaseSelector > > _reg_vtx_selectors
description
Definition VertexAnaProcessor.h:124
VertexAnaProcessor::reg_it
std::map< std::string, std::shared_ptr< TrackHistos > >::iterator reg_it
description
Definition VertexAnaProcessor.h:131
VertexAnaProcessor::vtxSelector
std::shared_ptr< BaseSelector > vtxSelector
description
Definition VertexAnaProcessor.h:88
VertexAnaProcessor::_vtx_histos
std::shared_ptr< TrackHistos > _vtx_histos
description
Definition VertexAnaProcessor.h:120
VertexAnaProcessor::ecal_
std::vector< CalCluster * > * ecal_
description
Definition VertexAnaProcessor.h:103
VertexAnaProcessor::evth_
EventHeader * evth_
description
Definition VertexAnaProcessor.h:101
VertexAnaProcessor::bmcParts_
TBranch * bmcParts_
description
Definition VertexAnaProcessor.h:97
VertexAnaProcessor::hits_
std::vector< TrackerHit * > * hits_
description
Definition VertexAnaProcessor.h:106
VertexAnaProcessor::tsColl_
std::string tsColl_
description
Definition VertexAnaProcessor.h:110
VertexAnaProcessor::v0proj_fits_
json v0proj_fits_
json object v0proj
Definition VertexAnaProcessor.h:148
VertexAnaProcessor::beamPosCfg_
std::string beamPosCfg_
json containing run dep beamspot positions
Definition VertexAnaProcessor.h:144
VertexAnaProcessor::_reg_tuples
std::map< std::string, std::shared_ptr< FlatTupleMaker > > _reg_tuples
description
Definition VertexAnaProcessor.h:127
VertexAnaProcessor::_ah
std::shared_ptr< AnaHelpers > _ah
description
Definition VertexAnaProcessor.h:141
VertexAnaProcessor::bevth_
TBranch * bevth_
description
Definition VertexAnaProcessor.h:98
VertexAnaProcessor::trkColl_
std::string trkColl_
description
Definition VertexAnaProcessor.h:113
VertexAnaProcessor::brMap_
std::map< const char *, int, char_cmp > brMap_
description
Definition VertexAnaProcessor.h:92
VertexAnaProcessor::btrks_
TBranch * btrks_
description
Definition VertexAnaProcessor.h:96
VertexAnaProcessor::vtxColl_
std::string vtxColl_
description
Definition VertexAnaProcessor.h:111
VertexAnaProcessor::ecalColl_
std::string ecalColl_
description
Definition VertexAnaProcessor.h:114
VertexAnaProcessor::_reg_mc_vtx_histos
std::map< std::string, std::shared_ptr< MCAnaHistos > > _reg_mc_vtx_histos
description
Definition VertexAnaProcessor.h:126
VertexAnaProcessor::bts_
TBranch * bts_
description
Definition VertexAnaProcessor.h:93
VertexAnaProcessor::hitColl_
std::string hitColl_
description
Definition VertexAnaProcessor.h:112
VertexAnaProcessor::mcHistoCfg_
std::string mcHistoCfg_
description
Definition VertexAnaProcessor.h:135
VertexAnaProcessor::vtxs_
std::vector< Vertex * > * vtxs_
description
Definition VertexAnaProcessor.h:104
VertexAnaProcessor::beamPosCorrections_
std::vector< double > beamPosCorrections_
holds beam position corrections
Definition VertexAnaProcessor.h:146
VertexAnaProcessor::bhits_
TBranch * bhits_
description
Definition VertexAnaProcessor.h:95
VertexAnaProcessor::_regions
std::vector< std::string > _regions
description
Definition VertexAnaProcessor.h:129
VertexAnaProcessor::mcParts_
std::vector< MCParticle * > * mcParts_
description
Definition VertexAnaProcessor.h:107
utils::getKalmanTrackL1Isolations
double getKalmanTrackL1Isolations(Track *track, std::vector< TrackerHit * > *siClusters)
description
Definition utilities.cxx:566
utils::v0_projection_to_target_significance
double v0_projection_to_target_significance(json v0proj_fits, int run, double &vtx_proj_x, double &vtx_proj_y, double &vtx_proj_x_signif, double &vtx_proj_y_signif, double vtx_x, double vtx_y, double vtx_z, double vtx_px, double vtx_py, double vtx_pz)
description
Definition utilities.cxx:879
utils::get2016KFMCTruthHitCodes
void get2016KFMCTruthHitCodes(Track *ele_trk, Track *pos_trk, int &L1L2hitCode, int &L1hitCode, int &L2hitCode)
description
Definition utilities.cxx:675
TSData::tsBits::Single_2_Bot
bool Single_2_Bot
Definition TSData.h:49
TSData::tsBits::Single_2_Top
bool Single_2_Top
Definition TSData.h:45
TSData::tsBits::Single_3_Top
bool Single_3_Top
Definition TSData.h:46
TSData::tsBits::Single_3_Bot
bool Single_3_Bot
Definition TSData.h:50
Generated by doxygen 1.9.8