Class Reference for E1039 Core & Analysis Software
AnaSortMixVertex.cc
Go to the documentation of this file.
1 #include <fstream>
2 #include <ctime>
3 #include <iomanip>
4 #include <TCanvas.h>
5 #include <TH1D.h>
6 #include <TSystem.h>
7 #include <TFile.h>
8 #include <TTree.h>
9 #include <phool/recoConsts.h>
10 #include <geom_svc/GeomSvc.h>
11 #include <ktracker/SRawEvent.h>
12 #include <ktracker/SRecEvent.h>
13 #include <interface_main/SQTrackVector.h>
14 #include <GenFit/FieldManager.h>
15 #include <phfield/PHFieldConfig.h>
16 #include <phfield/PHFieldConfig_v3.h>
17 #include <phfield/PHFieldUtility.h>
18 #include "AnaSortMixVertex.h"
19 #include "SQVertexing_v2.h"
20 #include <UtilAna/UtilDimuon.h>
21 #include "UtilBeam.h"
22 #include "TreeData.h"
23 
24 using namespace std;
25 
26 AnaSortMixVertex::AnaSortMixVertex(const std::string name)
27  : m_name(name)
28  , m_dir_out(".")
29  , m_ds("")
30  , m_run_id(0)
31  , m_use_raw(true)
32  , m_use_rec(true)
33  , m_n_evt_max(0)
34  , m_verb(0)
35  , m_raw(0)
36  , m_rec(0)
37  , ana_tree(0)
38  , mixed_tree(0)
39  , m_target_pos(-1)
40 {
41  ;
42 }
43 AnaSortMixVertex::~AnaSortMixVertex()
44 {
45  if (m_raw) delete m_raw;
46  if (m_rec) delete m_rec;
47 }
48 
50 
53 void AnaSortMixVertex::Init(const int run_id)
54 {
55  //initializing the reco constants used during original reconstruction
56  recoConsts* rc = recoConsts::instance();
57  rc->set_IntFlag("RUNNUMBER", run_id);
58  rc->set_DoubleFlag("FMAGSTR", -1.044);
59  rc->set_DoubleFlag("KMAGSTR", -1.025);
60  rc->set_DoubleFlag("X_BEAM", -0.2);
61  rc->set_BoolFlag("COARSE_MODE", false);
62  rc->set_BoolFlag("REQUIRE_MUID", false);
63  rc->set_CharFlag("HIT_MASK_MODE", "X");
64 
65  rc->set_CharFlag("AlignmentMille", "$E1039_RESOURCE/alignment/run0/align_mille_v10_a.txt");
66  rc->set_CharFlag("AlignmentHodo", "");
67  rc->set_CharFlag("AlignmentProp", "");
68  rc->set_CharFlag("Calibration", "");
69  rc->set_IntFlag ("MaxHitsDC0" , 350);
70  rc->set_IntFlag ("MaxHitsDC1" , 350);
71  rc->set_IntFlag ("MaxHitsDC2" , 170);
72  rc->set_IntFlag ("MaxHitsDC3p", 140);
73  rc->set_IntFlag ("MaxHitsDC3m", 140);
74  rc->set_DoubleFlag("RejectWinDC0" , 0.3);
75  rc->set_DoubleFlag("RejectWinDC1" , 0.5);
76  rc->set_DoubleFlag("RejectWinDC2" , 0.35);
77  rc->set_DoubleFlag("RejectWinDC3p", 0.24);
78  rc->set_DoubleFlag("RejectWinDC3m", 0.24);
79 
80  //setting up the output directory
81  if (m_dir_out != ".") gSystem->mkdir(m_dir_out.c_str(), true);
82  m_ofs.open((m_dir_out + "/output.txt").c_str());
83  m_file_out = new TFile((m_dir_out + "/output.root").c_str(), "RECREATE");
84 }
85 
86 void AnaSortMixVertex::Analyze(const int run_id, const std::vector<std::string> list_in)
87 {
88  //this is the main funtion to collect all the data from all good spill for given run and call the sort, mix and vertec function.
89 
90  //tree to store all the necessary information collected from the input reco files
91  m_file_out->cd();
92  ana_tree = new TTree("ana_tree","tree to store analysis variables");
93  ana_tree->Branch("pos_tracks", &pos_tracks);
94  ana_tree->Branch("neg_tracks", &neg_tracks);
95  ana_tree->Branch("fpga1", &fpga1, "fpga1/I");
96  ana_tree->Branch("fpga2", &fpga2, "fpga2/I");
97  ana_tree->Branch("fpga3", &fpga3, "fpga3/I");
98  ana_tree->Branch("fpga4", &fpga4, "fpga4/I");
99  //ana_tree->Branch("occuD0", &occuD0, "occuD0/I");
100  ana_tree->Branch("occuD1", &occuD1, "occuD1/I"); //using occuD1 to store occuD0 for e1039.
101  ana_tree->Branch("occuD2", &occuD2, "occuD2/I");
102  ana_tree->Branch("occuD3p", &occuD3p, "occuD3p/I");
103  ana_tree->Branch("occuD3m", &occuD3m, "occuD3m/I");
104  ana_tree->Branch("runID", &runID, "runID/I");
105  ana_tree->Branch("spill_ID", &spill_ID, "spill_ID/I");
106  ana_tree->Branch("event_ID", &event_ID, "event_ID/I");
107  ana_tree->Branch("TargetPos", &TargetPos, "TargetPos/I");
108  ana_tree->Branch("rfp00c", &rfp00c, "rfp00c/f");
109  ana_tree->Branch("pot_p00", &pot_p00, "pot_p00/f");
110  ana_tree->Branch("liveP", &liveP, "liveP/f");
111 
112  TFile* file_rec;
113  TTree* tree_rec;
114 
115  cout<<"\n\n================= Collecting track and event info starts ===================" <<endl;
116 
117  std::clock_t start;
118  double duration;
119  start = std::clock();
120 
121  int events_tot = 0;
122 
123  for (auto it = list_in.begin(); it != list_in.end(); it++) {
124  cout << "Input = " << *it << endl;
125  const string fn_raw = *it;
126  const string fn_rec = *it;
127 
128  unsigned int n_evt_rec = 0;
129 
130  file_rec = new TFile(fn_rec.c_str());
131  if (! file_rec) {
132  cout << "Cannot get the DST tree. Continue." << endl;
133  continue;
134  //exit(1);
135  }
136 
137  tree_rec = (TTree*)file_rec->Get("T");
138  if (! tree_rec) {
139  cout << "Cannot find the rec data tree. Continue." << endl;
140  continue;
141  //exit(1);
142  }
143  if (! m_rec) m_rec = new SRecEvent();
144 
145  SQEvent *sq_evt = 0;
146  SQHitVector *sq_hitvec = 0;
147  SQTrackVector* sq_trk_vec = 0;
148  SRawEvent *m_org = new SRawEvent();
149  tree_rec->SetBranchAddress("DST.SQEvent", &sq_evt);
150  tree_rec->SetBranchAddress("DST.SQRecTrackVector", &sq_trk_vec);
151  //tree_rec->SetBranchAddress("DST.SRecEvent", &m_rec);
152  tree_rec->SetBranchAddress("DST.SQHitVector", &sq_hitvec);
153  n_evt_rec = tree_rec->GetEntries();
154 
155  cout << "N of events: all = " << n_evt_rec << ", to be analyzed = " << n_evt_rec << endl;
156 
157 
158 
160 
161 
162  for (std::size_t i_evt = 0; i_evt < n_evt_rec; i_evt++) {
163 
164  if ((i_evt+1) % 10000 == 0) cout << i_evt+1 << endl;
165  else if ((i_evt+1) % 1000 == 0) cout << " . " << flush;
166  //tree_raw->GetEntry(i_evt);
167  tree_rec->GetEntry(i_evt);
168  fpga1=0;
169  fpga2=0;
170  fpga3=0;
171  fpga4=0;
172 
173  //occuD0=0; for now i will store D0 occu in occuD1 so I do not have to modify the rest of the code.
174  occuD1=0;
175  occuD2=0;
176  occuD3p=0;
177  occuD3m=0;
178  TargetPos=-9999;
179  spill_ID = -9999;
180  event_ID = 0;
181  rfp00c = 0.;
182  pot_p00 = 0.;
183  liveP = 0.; //Live PoT
184 
185  pos_tracks.clear();
186  neg_tracks.clear();
187 
188  UtilSRawEvent::SetEvent (m_org, sq_evt, true);
189  UtilSRawEvent::SetHit (m_org, sq_hitvec);
190 
191  runID = m_org->getRunID();
192  spill_ID = m_org->getSpillID();
193  event_ID = m_org->getEventID();
194 
195 
196 
197  liveP = UtilBeam::PoTLive(spill_ID);
198  TargetPos = m_org->getTargetPos();
199 
200  if(m_target_pos>0 && TargetPos != m_target_pos) continue;
201  if (!m_org->isTriggeredBy(SRawEvent::MATRIX1)) continue;
202  events_tot ++;
203  fpga1 = (m_org->isTriggeredBy(SRawEvent::MATRIX1));
204  fpga2 = (m_org->isTriggeredBy(SRawEvent::MATRIX2));
205  fpga3 = (m_org->isTriggeredBy(SRawEvent::MATRIX3));
206  fpga4 = (m_org->isTriggeredBy(SRawEvent::MATRIX4));
207 
208  int rfp00 = m_org->getIntensity(0);
209  //double qie_ped = e906sc->GetQIEPedestal(spill_ID);
210  double qie_ped = 0; //set to zero for now as we do not have a datacatalog for e1039
211  rfp00c = rfp00 - qie_ped;
212  double coef = UtilBeam::PoTPerQIE(spill_ID);
213  pot_p00 = coef * rfp00c;
214 
215  occuD1 = m_org->getNHitsInD0();
216  occuD2 = m_org->getNHitsInD2();
217  occuD3p = m_org->getNHitsInD3p();
218  occuD3m = m_org->getNHitsInD3m();
219 
220  for (auto it = sq_trk_vec->begin(); it != sq_trk_vec->end(); it++) {
221  //int TrkMult = m_rec->getNTracks();
222  //Loop to get vector<SRecTrack> for +ve and -ve tracks from rec event
223  //for (int i =0; i< TrkMult; i++)
224  //{
225  //SRecTrack recTrack = m_rec->getTrack(i);
226  SRecTrack* recTrack = dynamic_cast<SRecTrack*>(*it);
227  if (recTrack->getCharge()>0)
228  {
229  pos_tracks.push_back(*recTrack);
230  }
231  if (recTrack->getCharge()<0)
232  {
233  neg_tracks.push_back(*recTrack);
234  }
235 
236  }
237 
238  ana_tree->Fill();
239 
240 
241  }
242  std::cout<<std::endl;
243  if (m_use_rec) file_rec->Close();
244  std::cout<<"Current total events in the run = " <<events_tot<<std::endl;
245  }
246  m_file_out->cd();
247  ana_tree->Write();
248  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
249 
250  std::cout<<"\n================Collecting track and envent info Ended, cpu time taken: ==="<< duration <<endl;
251 
254  //DoVertex(ana_tree, run_id, 1);
255  SortTree(ana_tree); //sort the tree
256  MixTracks(sorted_tree);//Mix the sorted tree
257  DoVertex(sorted_tree, run_id, 0);//Vertex the sorted tree
258  DoVertex(mixed_tree, run_id, 1);//Vertex the mixed tree
259 
262 
263  cout<< endl;
264 }
265 
266 
267 void AnaSortMixVertex::End()
268 {
269  m_ofs.close();
270  m_file_out->cd();
271  m_file_out->Close();
272 }
273 
274 
275 
277 
281 void AnaSortMixVertex::SortTree(TTree* tree_to_sort)
282 {
283 
284  cout<<"\n\n================= Sorting starts ===================" <<endl;
285 
286  std::clock_t start;
287  double duration;
288  start = std::clock();
289 
290  m_file_out->cd();
291  sorted_tree = (TTree*)tree_to_sort->CloneTree(0);
292 
293 
294  Int_t nentries = (Int_t)tree_to_sort->GetEntries();
295 
296  tree_to_sort->Draw("occuD1","","goff");
297 
298  Int_t *index = new Int_t[nentries];
299  //sort array containing occuD1 in decreasing order
300  //The array index contains the entry numbers in decreasing order of occuD1
301  TMath::Sort(nentries,tree_to_sort->GetV1(),index);
302 
303 
304  cout<<" No. of enteries to be sorted: "<<nentries<<endl;
305  for (Int_t i=0;i<nentries;i++) {
306  tree_to_sort->GetEntry(index[i]);
307  tree_to_sort->LoadBaskets(2000000000);
308  sorted_tree->Fill();
309  }
310  tree_to_sort->DropBaskets();
311  //tree_to_sort = NULL; //pointed dangling ptr to NULL
312  sorted_tree->SetName("sorted_tree");
313  tree_to_sort->Delete();
314  m_file_out->cd();
315  sorted_tree->Write();
316 
317  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
318  std::cout<<"\n================ Sorting Ended, cpu time taken: ==="<< duration <<endl;
319 
320  delete [] index;
321 
322 }
323 
324 
326 
327 void AnaSortMixVertex::MixTracks(TTree* sorted_tree1)
328 {
329 
330  cout<<"\n\n====================== Mixing Begins========================"<<endl;
331  std::clock_t start;
332  double duration;
333  start = std::clock();
334 
335  //tree to store mixed tracks and their information
336  m_file_out->cd();
337  mixed_tree = new TTree("mixed_tree", "tree from mixing");
338  mixed_tree->Branch("pos_tracks", &pos_tracks_mix);
339  mixed_tree->Branch("neg_tracks", &neg_tracks_mix);
340 
341  mixed_tree->Branch("plus_fpga1", &plus_fpga1);
342  mixed_tree->Branch("minus_fpga1", &minus_fpga1);
343 
344  mixed_tree->Branch("plus_occuD1", &plus_occuD1, "plus_occuD1/I");
345  mixed_tree->Branch("plus_occuD2", &plus_occuD2, "plus_occuD2/I");
346  mixed_tree->Branch("plus_occuD3m", &plus_occuD3m, "plus_occuD3m/I");
347  mixed_tree->Branch("plus_occuD3p", &plus_occuD3p, "plus_occuD3p/I");
348  mixed_tree->Branch("plus_rfp00c", &plus_rfp00c, "plus_rfp00c/f");
349  mixed_tree->Branch("plus_pot_p00", &plus_pot_p00, "plus_pot_p00/f");
350  mixed_tree->Branch("plus_liveP", &plus_liveP, "plus_liveP/f");
351  mixed_tree->Branch("plus_TargetPos", &plus_TargetPos, "plus_TargetPos/I");
352  mixed_tree->Branch("plus_spill_ID", &plus_spill_ID, "plus_spill_ID/I");
353  mixed_tree->Branch("plus_event_ID", &plus_event_ID, "plus_event_ID/I");
354 
355 
356  mixed_tree->Branch("minus_occuD1", &minus_occuD1, "minus_occuD1/I");
357  mixed_tree->Branch("minus_occuD2", &minus_occuD2, "minus_occuD2/I");
358  mixed_tree->Branch("minus_occuD3m", &minus_occuD3m, "minus_occuD3m/I");
359  mixed_tree->Branch("minus_occuD3p", &minus_occuD3p, "minus_occuD3p/I");
360  mixed_tree->Branch("minus_rfp00c", &minus_rfp00c, "minus_rfp00c/f");
361  mixed_tree->Branch("minus_pot_p00", &minus_pot_p00, "minus_pot_p00/f");
362  mixed_tree->Branch("minus_liveP", &minus_liveP, "minus_liveP/f");
363  mixed_tree->Branch("minus_TargetPos", &minus_TargetPos, "minus_TargetPos/I");
364  mixed_tree->Branch("minus_spill_ID", &minus_spill_ID, "minus_spill_ID/I");
365  mixed_tree->Branch("minus_event_ID", &minus_event_ID, "minus_event_ID/I");
366 
367 
368  std::vector<SRecTrack>* pos_tracks_ = 0;
369  std::vector<SRecTrack>* neg_tracks_ = 0;
370 
371  int occuD1_,occuD2_, occuD3p_, occuD3m_,TargetPos_, fpga1_, event_ID_, spill_ID_;
372  occuD1_ = occuD2_ = occuD3p_ = occuD3m_ = TargetPos_ = event_ID_ = spill_ID_ = 0;
373  float rfp00c_, pot_p00_, liveP_;
374  rfp00c_ = pot_p00_ = liveP_ = 0.0;
375 
376  //getting the input tree branches
377  sorted_tree1->SetBranchAddress("pos_tracks", &pos_tracks_);
378  sorted_tree1->SetBranchAddress("neg_tracks", &neg_tracks_);
379  sorted_tree1->SetBranchAddress("fpga1", &fpga1_);
380  sorted_tree1->SetBranchAddress("occuD1", &occuD1_);
381  sorted_tree1->SetBranchAddress("occuD2", &occuD2_);
382  sorted_tree1->SetBranchAddress("occuD3p", &occuD3p_);
383  sorted_tree1->SetBranchAddress("occuD3m", &occuD3m_);
384  sorted_tree1->SetBranchAddress("TargetPos", &TargetPos_);
385  sorted_tree1->SetBranchAddress("rfp00c", &rfp00c_);
386  sorted_tree1->SetBranchAddress("pot_p00", &pot_p00_);
387  sorted_tree1->SetBranchAddress("liveP", &liveP_);
388  sorted_tree1->SetBranchAddress("spill_ID", &spill_ID_);
389  sorted_tree1->SetBranchAddress("event_ID", &event_ID_);
390 
391  cout<<"No. of entries from sorted tree: "<<sorted_tree1->GetEntries()<<endl;
392 
393  for(int i_evt=0;i_evt<sorted_tree1->GetEntries();i_evt++)
394  {
395 
399  sorted_tree1->GetEntry(i_evt);
400  plus_fpga1 = fpga1_;
401  plus_occuD1 = occuD1_;
402  plus_occuD2 = occuD2_;
403  plus_occuD3p = occuD3p_;
404  plus_occuD3m = occuD3m_;
405  plus_TargetPos = TargetPos_;
406  plus_rfp00c = rfp00c_;
407  plus_pot_p00 = pot_p00_;
408  plus_liveP = liveP_;
409  plus_event_ID = event_ID_;
410  plus_spill_ID = spill_ID_;
411 
412  for (std::size_t j=0; j<pos_tracks_->size(); ++j)
413  {
414  pos_tracks_mix.push_back(pos_tracks_->at(j));
415  }
416 
420  sorted_tree1->GetEntry(i_evt+1);
421  minus_fpga1 = fpga1_;
422  minus_occuD1 = occuD1_;
423  minus_occuD2 = occuD2_;
424  minus_occuD3p = occuD3p_;
425  minus_occuD3m = occuD3m_;
426  minus_TargetPos = TargetPos_;
427  minus_rfp00c = rfp00c_;
428  minus_pot_p00 = pot_p00_;
429  minus_liveP = liveP_;
430  minus_event_ID = event_ID_;
431  minus_spill_ID = spill_ID_;
432 
433  for (std::size_t k=0; k<neg_tracks_->size(); ++k)
434  {
435  neg_tracks_mix.push_back(neg_tracks_->at(k));
436  }
437 
438  mixed_tree->Fill();
439  pos_tracks_mix.clear();
440  neg_tracks_mix.clear();
441 
442  plus_occuD1 = 0;
443  plus_occuD2 = 0;
444  plus_occuD3p = 0;
445  plus_occuD3m = 0;
446  plus_TargetPos = 0;
447  plus_rfp00c = 0;
448  plus_pot_p00 = 0;
449  plus_liveP = 0;
450  plus_event_ID = 0;
451  plus_spill_ID = 0;
452 
453  minus_occuD1 = 0;
454  minus_occuD2 = 0;
455  minus_occuD3p = 0;
456  minus_occuD3m = 0;
457  minus_TargetPos = 0;
458  minus_rfp00c = 0;
459  minus_pot_p00 = 0;
460  minus_liveP = 0;
461  minus_event_ID = 0;
462  minus_spill_ID = 0;
463  }
464 
465  m_file_out->cd();
466  mixed_tree->Write();
467  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
468  std::cout<<"\n================ Mixing Ended, CPU time taken in sec: ==="<< duration<<endl;
469 
470 }
471 
475 void AnaSortMixVertex::DoVertex(TTree* input_tree, const int run_id, bool mix_flag)
476 {
477 
478  cout<<"\n\n======================== Vertexing Starts ================ "<<endl;
479  cout<<"Mixflag: "<<mix_flag<<endl;
480  std::clock_t start;
481  double duration;
482  start = std::clock();
483 
484  std::cout<<"No. of events in the tree "<<input_tree->GetEntries()<<std::endl;
485 
486  recEvent = new SRecEvent();
487 
488  m_file_out->cd();
489  if(!mix_flag)
490  {
491  //tree to store results after vertexing the sorted_tree
492  save_sorted = new TTree("save_sorted", "the sorted tree after vertexing");
493  save_sorted->Branch("recEvent", &recEvent, 256000, 99);
494  save_sorted->Branch("fpga1", &fpga1, "fpga1/I");
495  save_sorted->Branch("occuD1", &occuD1, "occuD1/I");
496  save_sorted->Branch("occuD2", &occuD2, "occuD2/I");
497  save_sorted->Branch("occuD3m", &occuD3m, "occuD3m/I");
498  save_sorted->Branch("occuD3p", &occuD3p, "occuD3p/I");
499  save_sorted->Branch("rfp00c", &rfp00c, "rfp00c/f");
500  save_sorted->Branch("pot_p00", &pot_p00, "pot_p00/f");
501  save_sorted->Branch("liveP", &liveP, "liveP/f");
502  save_sorted->Branch("TargetPos", &TargetPos, "TargetPos/I");
503  save_sorted->Branch("spill_ID", &spill_ID, "spill_ID/I");
504  save_sorted->Branch("event_ID", &event_ID, "event_ID/I");
505  }
506  else
507  {
508  //tree to store results after vertexing the mixed_tree
509  save_mix = new TTree("save_mix", "the mixed tree after vertexing");
510  save_mix->Branch("recEvent", &recEvent, 256000, 99);
511  save_mix->Branch("plus_fpga1", &plus_fpga1, "plus_fpga1/I");
512  save_mix->Branch("plus_occuD1", &plus_occuD1, "plus_occuD1/I");
513  save_mix->Branch("plus_occuD2", &plus_occuD2, "plus_occuD2/I");
514  save_mix->Branch("plus_occuD3m", &plus_occuD3m, "plus_occuD3m/I");
515  save_mix->Branch("plus_occuD3p", &plus_occuD3p, "plus_occuD3p/I");
516  save_mix->Branch("plus_rfp00c", &plus_rfp00c, "plus_rfp00c/f");
517  save_mix->Branch("plus_pot_p00", &plus_pot_p00, "plus_pot_p00/f");
518  save_mix->Branch("plus_liveP", &plus_liveP, "plus_liveP/f");
519  save_mix->Branch("plus_TargetPos", &plus_TargetPos, "plus_TargetPos/I");
520  save_mix->Branch("plus_spill_ID", &plus_spill_ID, "plus_spill_ID/I");
521  save_mix->Branch("plus_event_ID", &plus_event_ID, "plus_event_ID/I");
522 
523  save_mix->Branch("minus_fpga1", &minus_fpga1, "minus_fpga1/I");
524  save_mix->Branch("minus_occuD1", &minus_occuD1, "minus_occuD1/I");
525  save_mix->Branch("minus_occuD2", &minus_occuD2, "minus_occuD2/I");
526  save_mix->Branch("minus_occuD3m", &minus_occuD3m, "minus_occuD3m/I");
527  save_mix->Branch("minus_occuD3p", &minus_occuD3p, "minus_occuD3p/I");
528  save_mix->Branch("minus_rfp00c", &minus_rfp00c, "minus_rfp00c/f");
529  save_mix->Branch("minus_pot_p00", &minus_pot_p00, "minus_pot_p00/f");
530  save_mix->Branch("minus_liveP", &minus_liveP, "minus_liveP/f");
531  save_mix->Branch("minus_TargetPos", &minus_TargetPos, "minus_TargetPos/I");
532  save_mix->Branch("minus_spill_ID", &minus_spill_ID, "minus_spill_ID/I");
533  save_mix->Branch("minus_event_ID", &minus_event_ID, "minus_event_ID/I");
534 
535  }
536 
537  //VertexFit* vtxfit = new VertexFit();
538  SQVertexing_v2* vtxfit = new SQVertexing_v2();
539  vtxfit->set_geom_file_name((string)gSystem->Getenv("E1039_RESOURCE") + "/geometry/geom_run005433.root");
540  vtxfit->set_legacy_rec_container(true);
541  gfield_v2 = nullptr;
542  vtxfit->InitVar();
543  //vtxfit->Verbosity(99);
544  vtxfit->InitField();
545  vtxfit->InitGeom();
546 
547  SRecDimuon recDimuon;
548 
549  std::vector<SRecTrack> * pos_tracks_input = 0;
550  std::vector<SRecTrack> * neg_tracks_input = 0;
551 
552  int plus_fpga1_, plus_occuD1_,plus_occuD2_,plus_occuD3p_, plus_occuD3m_, plus_TargetPos_;
553  plus_occuD1_ =plus_occuD2_ = plus_occuD3p_ = plus_occuD3m_ = plus_TargetPos_ = 0;
554  float plus_rfp00c_, plus_pot_p00_, plus_liveP_ ;
555  plus_rfp00c_ = plus_pot_p00_ = plus_liveP_ = 0.0;
556 
557  int minus_fpga1_, minus_occuD1_, minus_occuD2_, minus_occuD3p_,minus_occuD3m_,minus_TargetPos_;
558  minus_occuD1_ = minus_occuD2_ = minus_occuD3p_ = minus_occuD3m_ = minus_TargetPos_ = 0;
559  float minus_rfp00c_, minus_pot_p00_, minus_liveP_;
560  minus_rfp00c_ = minus_pot_p00_ = minus_liveP_ = 0.;
561 
562  int fpga1_, occuD1_, occuD2_, occuD3p_, occuD3m_, TargetPos_;
563  int event_ID_, spill_ID_, plus_event_ID_, plus_spill_ID_,minus_event_ID_, minus_spill_ID_;
564  occuD1_ = occuD2_ = occuD3p_ = occuD3m_ =TargetPos_ = 0;
565  event_ID_=spill_ID_=plus_event_ID_ = plus_spill_ID_ = minus_event_ID_ = minus_spill_ID_= 0;
566 
567  float rfp00c_, pot_p00_ , liveP_;
568  rfp00c_ = pot_p00_ = liveP_ = 0;
569 
570  input_tree->SetBranchAddress("pos_tracks", &pos_tracks_input);
571  input_tree->SetBranchAddress("neg_tracks", &neg_tracks_input);
572 
573  if(mix_flag)
574  {
575  input_tree->SetBranchAddress("plus_fpga1", &plus_fpga1_);
576  input_tree->SetBranchAddress("plus_occuD1", &plus_occuD1_);
577  input_tree->SetBranchAddress("plus_occuD2", &plus_occuD2_);
578  input_tree->SetBranchAddress("plus_occuD3m", &plus_occuD3m_);
579  input_tree->SetBranchAddress("plus_occuD3p", &plus_occuD3p_);
580  input_tree->SetBranchAddress("plus_rfp00c", &plus_rfp00c_);
581  input_tree->SetBranchAddress("plus_pot_p00", &plus_pot_p00_);
582  input_tree->SetBranchAddress("plus_liveP", &plus_liveP_);
583  input_tree->SetBranchAddress("plus_TargetPos", &plus_TargetPos_);
584  input_tree->SetBranchAddress("plus_spill_ID", &plus_spill_ID_);
585  input_tree->SetBranchAddress("plus_event_ID", &plus_event_ID_);
586 
587  input_tree->SetBranchAddress("minus_fpga1", &minus_fpga1_);
588  input_tree->SetBranchAddress("minus_occuD1", &minus_occuD1_);
589  input_tree->SetBranchAddress("minus_occuD2", &minus_occuD2_);
590  input_tree->SetBranchAddress("minus_occuD3m", &minus_occuD3m_);
591  input_tree->SetBranchAddress("minus_occuD3p", &minus_occuD3p_);
592  input_tree->SetBranchAddress("minus_rfp00c", &minus_rfp00c_);
593  input_tree->SetBranchAddress("minus_pot_p00", &minus_pot_p00_);
594  input_tree->SetBranchAddress("minus_liveP", &minus_liveP_);
595  input_tree->SetBranchAddress("minus_TargetPos", &minus_TargetPos_);
596  input_tree->SetBranchAddress("minus_spill_ID", &minus_spill_ID_);
597  input_tree->SetBranchAddress("minus_event_ID", &minus_event_ID_);
598  }
599  else
600  {
601  input_tree->SetBranchAddress("fpga1", &fpga1_);
602  input_tree->SetBranchAddress("occuD1", &occuD1_);
603  input_tree->SetBranchAddress("occuD2", &occuD2_);
604  input_tree->SetBranchAddress("occuD3p", &occuD3p_);
605  input_tree->SetBranchAddress("occuD3m", &occuD3m_);
606  input_tree->SetBranchAddress("TargetPos", &TargetPos_);
607  input_tree->SetBranchAddress("rfp00c", &rfp00c_);
608  input_tree->SetBranchAddress("pot_p00", &pot_p00_);
609  input_tree->SetBranchAddress("liveP", &liveP_);
610  input_tree->SetBranchAddress("spill_ID", &spill_ID_);
611  input_tree->SetBranchAddress("event_ID", &event_ID_);
612  }
613 
614  for(int i_evt=0;i_evt<input_tree->GetEntries();i_evt++)
615  {
616  input_tree->GetEntry(i_evt);
617 
618  occuD1 = 0;
619 
620  plus_occuD1 = 0;
621  plus_occuD2 = 0;
622  plus_occuD3p = 0;
623  plus_occuD3m = 0;
624  plus_TargetPos = -999;
625  plus_rfp00c = 0;
626  plus_pot_p00 = 0;
627  plus_liveP = 0;
628  plus_spill_ID = 0;
629  plus_event_ID = 0;
630 
631  minus_occuD1 = 0;
632  minus_occuD2 = 0;
633  minus_occuD3p = 0;
634  minus_occuD3m = 0;
635  minus_TargetPos = -999;
636  minus_rfp00c = 0;
637  minus_pot_p00 = 0;
638  minus_liveP = 0;
639  minus_spill_ID = 0;
640  minus_event_ID = 0;
641 
642  spill_ID = 0;
643  event_ID = 0;
644 
645  if(!mix_flag) {
646  spill_ID = spill_ID_;
647  event_ID = event_ID_;
648  }
649 
650  plus_fpga1 = plus_fpga1_;
651  plus_occuD1 = plus_occuD1_;
652  plus_occuD2 = plus_occuD2_;
653  plus_occuD3p = plus_occuD3p_;
654  plus_occuD3m = plus_occuD3m_;
655  plus_TargetPos = plus_TargetPos_;
656  plus_rfp00c = plus_rfp00c_;
657  plus_pot_p00 = plus_pot_p00_;
658  plus_liveP = plus_liveP_;
659  plus_spill_ID = plus_spill_ID_;
660  plus_event_ID = plus_event_ID_;
661 
662  minus_fpga1 = plus_fpga1_;
663  minus_occuD1 = minus_occuD1_;
664  minus_occuD2 = minus_occuD2_;
665  minus_occuD3p = minus_occuD3p_;
666  minus_occuD3m = minus_occuD3m_;
667  minus_TargetPos = minus_TargetPos_;
668  minus_rfp00c = minus_rfp00c_;
669  minus_pot_p00 = minus_pot_p00_;
670  minus_liveP = minus_liveP_;
671  minus_spill_ID = minus_spill_ID_;
672  minus_event_ID = minus_event_ID_;
673 
674  occuD1 = 0;
675  occuD2 = 0;
676  occuD3p = 0;
677  occuD3m = 0;
678  TargetPos = -9999;
679  rfp00c = 0;
680  pot_p00 = 0;
681  liveP = 0;
682 
683  fpga1 = fpga1_;
684  occuD1 = occuD1_;
685  occuD2 = occuD2_;
686  occuD3p = occuD3p_;
687  occuD3m = occuD3m_;
688  TargetPos = TargetPos_;
689  rfp00c = rfp00c_;
690  pot_p00 = pot_p00_;
691  liveP = liveP_;
692 
693  if(i_evt%1000==0) cout<<"No. of events analyzed: "<<i_evt<<endl;
694 
695  std::vector<SRecTrack> pos_tracks_hold, neg_tracks_hold ;
696 
697  for (std::size_t j=0; j<pos_tracks_input->size(); ++j)
698  {
699  SRecTrack* trk_p = &(pos_tracks_input->at(j));
700  pos_tracks_hold.push_back(*trk_p);
701  recEvent->insertTrack(*trk_p);
702  }
703 
704 
705  for (std::size_t k=0; k<neg_tracks_input->size(); ++k)
706  {
707  SRecTrack* trk_n = &(neg_tracks_input->at(k));
708  neg_tracks_hold.push_back(*trk_n);
709  recEvent->insertTrack(*trk_n);
710  }
711 
712  for(std::size_t j=0; j<pos_tracks_hold.size(); ++j)
713  {
714  for (std::size_t k=0; k<neg_tracks_hold.size(); ++k)
715  {
716 
717 
718  //Implementing FPGA1 trigger roadset requirement
719  if((pos_tracks_hold.at(j).getTriggerRoad() > 0 && neg_tracks_hold.at(k).getTriggerRoad() > 0) || (pos_tracks_hold.at(j).getTriggerRoad() < 0 && neg_tracks_hold.at(k).getTriggerRoad() < 0)) continue;
720 
721 
722 
723  SRecTrack* track1 = &pos_tracks_hold.at(j);
724  SRecTrack* track2 = &neg_tracks_hold.at(k);
725  if(!track1->isKalmanFitted() || !track2->isKalmanFitted()) continue;
726 
727  if(vtxfit->processOneDimuon(track1, track2, recDimuon)){
728  if (track1->getStateVector(0).GetNcols()==0) std::cout<<"track1 Ncols ==0 "<<std::endl;
729  if (track2->getStateVector(0).GetNcols()==0) std::cout<<"track2 Ncols ==0 "<<std::endl;
730 
731  recDimuon.set_track_id_pos(static_cast<int>(j));
732  recDimuon.set_track_id_neg(static_cast<int>(k) + static_cast<int>(pos_tracks_hold.size()));
733  recDimuon.calcVariables();
734 
735  recEvent->insertDimuon(recDimuon);
736  }
737  }
738  }
739 
740  if(mix_flag) save_mix->Fill();
741  else save_sorted->Fill();
742  recEvent->clear();
743  pos_tracks_hold.clear();
744  neg_tracks_hold.clear();
745  }
746 
747  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
748  std::cout<<"\n================ Vertexing Ended, CPU time taken in sec: ==="<< duration<<endl;
749  m_file_out->cd();
750  if (mix_flag) save_mix->Write();
751  else save_sorted->Write();
752 }
753 
754 
755 
756 //currently no embedding is done in e1039. So, commenting out all embedding stuffs.
757 
759 /*
760 void AnaSortMixVertex::PrepareMC(char* fn_mctree)
761 {
762 
763 
765 
766  std::vector<SRecTrack> * pos_tracks = 0;
767  std::vector<SRecTrack> * neg_tracks = 0;
768  std::vector<double> dim_M_mc;
769 
770  TTree* mc_tree = new TTree("mc_tree", "Tree with MC for mixing");
771  mc_tree->Branch("pos_tracks",&pos_tracks);
772  mc_tree->Branch("neg_tracks",&neg_tracks);
773  mc_tree->Branch("dim_M_mc", &dim_M_mc);
774 
775 
777  SRecEvent* recEvent = new SRecEvent();
778  SRawMCEvent* raw = new SRawMCEvent();
779 
780  cout<<"Reading MC tree: "<<fn_mctree<<endl;
781 
782  TFile* mc_file = new TFile (fn_mctree, "read");
783  TTree* save = (TTree*)mc_file->Get("save");
784 
785  save->SetBranchAddress("recEvent",&recEvent);
786  save->SetBranchAddress("rawEvent",&raw);
787 
788 
789  for(int i_evt=0;i_evt<save->GetEntries();i_evt++)
790  {
791  save->GetEntry(i_evt);
792 
793  pos_tracks->clear();
794  neg_tracks->clear();
795  dim_M_mc.clear();
796 
797  if(raw->isTriggeredBy(SRawEvent::MATRIX1)) continue;
798 
799  if(recEvent->getNDimuons()<1) continue;
800  for (int i =0; i<recEvent->getNTracks(); i++)
801  {
802  SRecTrack recTrack = recEvent->getTrack(i);
803  if (recTrack.getCharge()>0) pos_tracks->push_back(recTrack);
804  if (recTrack.getCharge()<0) neg_tracks->push_back(recTrack);
805  }
806 
807  if((pos_tracks->at(0).getTriggerRoad() > 0 && neg_tracks->at(0).getTriggerRoad() > 0) || (pos_tracks->at(0).getTriggerRoad() < 0 && neg_tracks->at(0).getTriggerRoad() < 0)) continue;
808 
809 
810  SRecDimuon recDimuon_mc = recEvent->getDimuon(0);
811  dim_M_mc.push_back(recDimuon_mc.mass);
812 
813  mc_tree ->Fill();
814 
815  }
816 
817 
818 }
819 
820 
822 */
823 void AnaSortMixVertex::EmbedMCSignal( char *file_mixed, char* file_mcsignal)
824 {
825 
827  cout<<"creating tree to store embedded event"<<endl;
828  TTree* mixed_embed_tree = new TTree("tree","MC signals embeded on mixed tree");
829 
830  std::vector<SRecTrack> * pos_tracks = 0;
831  std::vector<SRecTrack> * neg_tracks = 0;
832  std::vector<double> dim_M_mixed, dim_M_mc;
833  mixed_embed_tree->Branch("pos_tracks",&pos_tracks);
834  mixed_embed_tree->Branch("neg_tracks",&neg_tracks);
835  //mixed_embed_tree->Branch("dim_M_mixed", &dim_M_mixed);
836  mixed_embed_tree->Branch("dim_M_mc", &dim_M_mc);
837 
839 
840  cout<<"Reading MC tree: "<<file_mcsignal<<endl;
841 
842 
843  std::vector<SRecTrack> * pos_tracks_mc = 0;
844  std::vector<SRecTrack> * neg_tracks_mc = 0;
845  std::vector<double> * dim_M_mc_ = 0;
846 
847  TFile* mc_file = new TFile (file_mcsignal, "read");
848  TTree* mc_tree = (TTree*)mc_file->Get("mc_tree");
849 
850 
851  mc_tree->SetBranchAddress("pos_tracks", &pos_tracks_mc);
852  mc_tree->SetBranchAddress("neg_tracks", &neg_tracks_mc);
853  mc_tree->SetBranchAddress("dim_M_mc", &dim_M_mc_);
854 
855 
856  //Reading "T" tree form fn_mixedtree
857  std::vector<SRecTrack> * pos_tracks_ = 0;
858  std::vector<SRecTrack> * neg_tracks_ = 0;
859  cout<<"Reading Mixed tree: "<<file_mixed<<endl;
860 
861  TFile* _mixedFile = new TFile (file_mixed,"read");
862  TTree* T = (TTree*)_mixedFile->Get("tree");
863 
864  T->SetBranchAddress("pos_tracks", &pos_tracks_);
865  T->SetBranchAddress("neg_tracks", &neg_tracks_);
866 
867 
868 
869  int i_mc_evt=0;
870  cout<<"Total entries in the mixed tree: "<<T->GetEntries()<<endl;
871 
872  for(int i_evt=0;i_evt<T->GetEntries();i_evt++)
873  {
874  T->GetEntry(i_evt);
875  pos_tracks->clear();
876  neg_tracks->clear();
877  dim_M_mc.clear();
878 
879 
881 
882  for (std::size_t j=0; j<pos_tracks_->size(); ++j)
883  {
884  pos_tracks->push_back(pos_tracks_->at(j));
885  }
886 
887 
890  for (std::size_t k=0; k<neg_tracks_->size(); ++k)
891  {
892  neg_tracks->push_back(neg_tracks_->at(k));
893  }
894 
895 
897 
898  if(i_evt%25==0)//embed every 30 mixed events
899  {
900  mc_tree->GetEntry(i_mc_evt);
901  pos_tracks->push_back(pos_tracks_mc->at(0));
902  neg_tracks->push_back(neg_tracks_mc->at(0));
903  dim_M_mc.push_back(dim_M_mc_->at(0));
904  //cout<<"mc_dimuon\t"<<dim_M_mc_->at(0)<<endl;
905  i_mc_evt++;
906 
907  }
908 
909 
910  mixed_embed_tree ->Fill();
911 
912 
913  }
914 
915 
916 }
917 
SQVertexing_v2.h
AnaSortMixVertex::liveP
float liveP
proton number corrsponding to RF+00 (pedestal corrected)
Definition: AnaSortMixVertex.h:126
AnaSortMixVertex::pos_tracks_mix
std::vector< SRecTrack > pos_tracks_mix
Definition: AnaSortMixVertex.h:135
AnaSortMixVertex::AnaSortMixVertex
AnaSortMixVertex(const std::string name="ana_mixvertex")
Definition: AnaSortMixVertex.cc:26
AnaSortMixVertex::End
virtual void End()
Definition: AnaSortMixVertex.cc:267
AnaSortMixVertex::Init
virtual void Init(const int run_id)
Function to initialize all variables for the 1st analysis step.
Definition: AnaSortMixVertex.cc:53
AnaSortMixVertex::DoVertex
virtual void DoVertex(TTree *inputtree, const int run_id, bool mix)
Definition: AnaSortMixVertex.cc:475
AnaSortMixVertex::neg_tracks_mix
std::vector< SRecTrack > neg_tracks_mix
Definition: AnaSortMixVertex.h:135
AnaSortMixVertex::gfield_v2
SQGenFit::GFField * gfield_v2
Definition: AnaSortMixVertex.h:140
AnaSortMixVertex::MixTracks
virtual void MixTracks(TTree *sorted_tree)
Mixing tracks from sorted_tree.
Definition: AnaSortMixVertex.cc:327
AnaSortMixVertex::pos_tracks
std::vector< SRecTrack > pos_tracks
Definition: AnaSortMixVertex.h:134
AnaSortMixVertex::~AnaSortMixVertex
virtual ~AnaSortMixVertex()
Definition: AnaSortMixVertex.cc:43
AnaSortMixVertex::Analyze
virtual void Analyze(const int run_id, const std::vector< std::string > list_in)
Definition: AnaSortMixVertex.cc:86
AnaSortMixVertex::rfp00c
float rfp00c
RF+00 corrected for the pedestal.
Definition: AnaSortMixVertex.h:124
AnaSortMixVertex::recEvent
SRecEvent * recEvent
Definition: AnaSortMixVertex.h:136
AnaSortMixVertex::neg_tracks
std::vector< SRecTrack > neg_tracks
Definition: AnaSortMixVertex.h:134
AnaSortMixVertex::SortTree
virtual void SortTree(TTree *tree_to_sort)
Definition: AnaSortMixVertex.cc:281
AnaSortMixVertex::EmbedMCSignal
virtual void EmbedMCSignal(char *file_mixed, char *file_mcsignal)
MC preperation for embedding.
Definition: AnaSortMixVertex.cc:823
PHFlag::set_IntFlag
virtual void set_IntFlag(const std::string &name, const int flag)
Definition: PHFlag.cc:145
PHFlag::set_BoolFlag
virtual void set_BoolFlag(const std::string &name, const bool flag)
Definition: PHFlag.cc:179
PHFlag::set_DoubleFlag
virtual void set_DoubleFlag(const std::string &name, const double flag)
Definition: PHFlag.cc:77
SQEvent
An SQ interface class to hold one event header.
Definition: SQEvent.h:17
SQHitVector
An SQ interface class to hold a list of SQHit objects.
Definition: SQHitVector.h:32
SQTrackVector
An SQ interface class to hold a list of SQTrack objects.
Definition: SQTrackVector.h:19
SQTrackVector::begin
virtual ConstIter begin() const =0
SQTrackVector::end
virtual ConstIter end() const =0
SQVertexing_v2
Definition: SQVertexing_v2.h:21
SQVertexing_v2::set_legacy_rec_container
void set_legacy_rec_container(const bool enable=true)
Definition: SQVertexing_v2.h:31
SQVertexing_v2::InitGeom
int InitGeom(PHCompositeNode *topNode)
SQVertexing_v2::set_geom_file_name
void set_geom_file_name(const std::string &geomFileName)
Definition: SQVertexing_v2.h:36
SQVertexing_v2::InitField
int InitField(PHCompositeNode *topNode)
SQVertexing_v2::processOneDimuon
bool processOneDimuon(SRecTrack *track1, SRecTrack *track2, SRecDimuon &dimuon)
Definition: SQVertexing_v2.cc:299
SQVertexing_v2::InitVar
void InitVar()
Definition: SQVertexing_v2.cc:80
SRawEvent::getNHitsInD2
Int_t getNHitsInD2()
Definition: SRawEvent.cxx:455
SRawEvent::getSpillID
Int_t getSpillID()
Definition: SRawEvent.h:151
SRawEvent::getEventID
Int_t getEventID()
Definition: SRawEvent.h:150
SRawEvent::getNHitsInD3p
Int_t getNHitsInD3p()
Definition: SRawEvent.cxx:471
SRawEvent::getNHitsInD3m
Int_t getNHitsInD3m()
Definition: SRawEvent.cxx:482
SRawEvent::getIntensity
Int_t getIntensity()
Definition: SRawEvent.h:201
SRawEvent::getNHitsInD0
Int_t getNHitsInD0()
Definition: SRawEvent.cxx:433
SRawEvent::getTargetPos
Int_t getTargetPos()
Definition: SRawEvent.h:195
SRawEvent::getRunID
Int_t getRunID()
Definition: SRawEvent.h:149
SRawEvent::isTriggeredBy
bool isTriggeredBy(Int_t trigger)
Definition: SRawEvent.h:177
SRecDimuon::calcVariables
void calcVariables(int choice=0)
Calculate the kinematic vairables, 0 = vertex, 1 = target, 2 = dump.
Definition: SRecEvent.cxx:633
SRecDimuon::set_track_id_pos
virtual void set_track_id_pos(const int a)
Definition: SRecEvent.h:321
SRecDimuon::set_track_id_neg
virtual void set_track_id_neg(const int a)
Definition: SRecEvent.h:324
SRecEvent::insertTrack
void insertTrack(SRecTrack trk)
Insert tracks.
Definition: SRecEvent.h:466
SRecEvent::insertDimuon
void insertDimuon(SRecDimuon dimuon)
Insert dimuon.
Definition: SRecEvent.h:470
SRecEvent::clear
void clear()
Clear everything.
Definition: SRecEvent.cxx:802
SRecTrack::getCharge
Int_t getCharge() const
Gets.
Definition: SRecEvent.h:101
SRecTrack::getStateVector
TMatrixD getStateVector(Int_t i)
Definition: SRecEvent.h:109
SRecTrack::isKalmanFitted
Bool_t isKalmanFitted()
Fit status.
Definition: SRecEvent.h:147
recoConsts::instance
static recoConsts * instance()
Definition: recoConsts.cc:7
recoConsts::set_CharFlag
virtual void set_CharFlag(const std::string &name, const std::string &flag)
overide the virtual function to expand the environmental variables
Definition: recoConsts.cc:21
UtilBeam::PoTPerQIE
double PoTPerQIE(const unsigned int spill_id)
Return PoT/QIEsum.
Definition: UtilBeam.cc:38
UtilBeam::PoTLive
double PoTLive(const unsigned int spill_id)
Return the live number of PoT.
Definition: UtilBeam.cc:19
UtilSRawEvent::SetEvent
bool SetEvent(SRawEvent *sraw, const SQEvent *evt, const bool do_assert=false)
Definition: UtilSRawEvent.cxx:10
UtilSRawEvent::SetHit
bool SetHit(SRawEvent *sraw, const SQHitVector *hit_vec, std::map< int, size_t > *hitID_idx=0, const bool do_assert=false)
Definition: UtilSRawEvent.cxx:57