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