Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PairMaker.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PairMaker.cc
1 
2 #include "PairMaker.h"
3 
4 #include <cstdlib>
5 #include <cstdio>
6 #include <iostream>
7 #include <iomanip>
8 #include <fstream>
9 
10 #include "TFile.h"
11 #include "TLorentzVector.h"
12 #include "TRandom2.h"
13 
14 #include <trackbase_historic/SvtxTrack.h>
15 #include <trackbase_historic/SvtxTrackMap.h>
16 #include <trackbase_historic/SvtxVertex.h>
17 #include <trackbase_historic/SvtxVertexMap.h>
18 #include <trackbase/TrkrDefs.h>
19 
20 #include <globalvertex/GlobalVertexMap.h>
21 #include <globalvertex/GlobalVertex.h>
22 
23 #include <calobase/RawClusterContainer.h>
24 #include <calobase/RawCluster.h>
25 #include <calobase/RawClusterv1.h>
26 
27 #include <phool/getClass.h>
28 #include <phool/recoConsts.h>
29 #include <phool/PHCompositeNode.h>
30 #include <phool/PHIODataNode.h>
31 #include <phool/PHNodeIterator.h>
32 #include <phool/PHRandomSeed.h>
33 #include <fun4all/Fun4AllReturnCodes.h>
34 
35 #include "sPHElectron.h"
36 #include "sPHElectronv1.h"
37 #include "sPHElectronPair.h"
38 #include "sPHElectronPairv1.h"
39 #include "sPHElectronPairContainer.h"
40 #include "sPHElectronPairContainerv1.h"
41 
42 #include "trackpidassoc/TrackPidAssoc.h"
43 
44 //#include <gsl/gsl_rng.h>
45 
46 typedef PHIODataNode<PHObject> PHObjectNode_t;
47 
48 using namespace std;
49 
50 //==============================================================
51 
52 PairMaker::PairMaker(const std::string &name, const std::string &filename) : SubsysReco(name)
53 {
54  _ZMIN = -15.;
55  _ZMAX = 15.;
56  _multbins.push_back(0.);
57 // _multbins.push_back(75.);
58 // _multbins.push_back(250.);
59 // _multbins.push_back(900.);
60  _multbins.push_back(3000.);
61  _multbins.push_back(9999.);
62  _min_buffer_depth = 10;
63  _max_buffer_depth = 50;
64  outnodename = "ElectronPairs";
65  _rng = nullptr;
66  EventNumber=0;
67 }
68 
69 //==============================================================
70 
71 int PairMaker::Init(PHCompositeNode *topNode)
72 {
73 
74  _rng = new TRandom2();
75  _rng->SetSeed(0);
76 
77  PHNodeIterator iter(topNode);
78  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
79  if (!dstNode)
80  {
81  cerr << PHWHERE << " ERROR: Can not find DST node." << endl;
82  return Fun4AllReturnCodes::ABORTEVENT;
83  }
84 
85  sPHElectronPairContainerv1* eePairs = new sPHElectronPairContainerv1();
86 
87  PHCompositeNode *PairNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", outnodename));
88  if (!PairNode)
89  {
90  PHObjectNode_t *PairNode = new PHIODataNode<PHObject>(eePairs,outnodename.c_str(),"PHObject");
91  dstNode->addNode(PairNode);
92  cout << PHWHERE << " INFO: added " << outnodename << endl;
93  }
94  else { cout << PHWHERE << " INFO: " << outnodename << " already exists." << endl; }
95 
96  std::cout << "PairMaker::Init ended." << endl;
97  return Fun4AllReturnCodes::EVENT_OK;
98 }
99 
100 //==============================================================
101 
102 int PairMaker::InitRun(PHCompositeNode *topNode)
103 {
104  return Fun4AllReturnCodes::EVENT_OK;
105 }
106 
107 //==============================================================
108 
109 int PairMaker::process_event(PHCompositeNode *topNode)
110 {
111  return process_event_test(topNode);
112 }
113 
114 //======================================================================
115 
116 int PairMaker::process_event_test(PHCompositeNode *topNode) {
117  EventNumber++;
118 
119 // int howoften = 1;
120 // if((EventNumber-1)%howoften==0) {
121  cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
122 // }
123 
124  if(EventNumber==1) topNode->print();
125 
126  sPHElectronPairContainerv1 *eePairs = findNode::getClass<sPHElectronPairContainerv1>(topNode,outnodename.c_str());
127  if(!eePairs) { cerr << outnodename << " not found!" << endl; }
128  else { cout << "Found " << outnodename << " node." << endl; }
129 
130  GlobalVertexMap *global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
131  if(!global_vtxmap) {
132  cerr << PHWHERE << " ERROR: Can not find GlobalVertexMap node." << endl;
133  return Fun4AllReturnCodes::ABORTEVENT;
134  }
135  //cout << "Number of GlobalVertexMap entries = " << global_vtxmap->size() << endl;
136 
137  SvtxVertexMap *vtxmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
138  if(!vtxmap) {
139  cout << "SvtxVertexMap node not found!" << endl;
140  return Fun4AllReturnCodes::ABORTEVENT;
141  }
142  //cout << "Number of SvtxVertexMap entries = " << vtxmap->size() << endl;
143 
144  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
145  if(!trackmap) {
146  cerr << PHWHERE << " ERROR: Can not find SvtxTrackMap node." << endl;
147  return Fun4AllReturnCodes::ABORTEVENT;
148  }
149 
150  RawClusterContainer* cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
151  if(!cemc_clusters) {
152  cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
153  return Fun4AllReturnCodes::ABORTEVENT;
154  }
155  else { cout << "FOUND CLUSTER_CEMC node." << endl; }
156 
157  int mycount = 0;
158  RawClusterContainer::ConstRange cluster_range = cemc_clusters->getClusters();
159  for (RawClusterContainer::ConstIterator cluster_iter = cluster_range.first; cluster_iter != cluster_range.second; cluster_iter++)
160  {
161  //RawCluster *cluster = cluster_iter->second;
162  //double phi = cluster->get_phi();
163  //double z = cluster->get_z();
164  //double ee = cluster->get_energy();
165  //int ntowers = cluster->getNTowers();
166  //if(ee>2.) cout << "cluster: " << ee << " " << ntowers << " " << phi << " " << z << endl;
167  mycount++;
168  }
169  cout << "Number of CEMC clusters = " << mycount << endl;
170 
171  TrackPidAssoc *track_pid_assoc = findNode::getClass<TrackPidAssoc>(topNode, "TrackPidAssoc");
172  if(!track_pid_assoc) {
173  cerr << PHWHERE << "ERROR: CAN NOT FIND TrackPidAssoc Node!" << endl;
174  return Fun4AllReturnCodes::ABORTEVENT;
175  }
176  auto electrons = track_pid_assoc->getTracks(TrackPidAssoc::electron);
177 
178 // for(auto it = electrons.first; it != electrons.second; ++it)
179 // {
180 // SvtxTrack *tr = trackmap->get(it->second);
181 // double p = tr->get_p();
182 // std::cout << " pid " << it->first << " track ID " << it->second << " mom " << p << std::endl;
183 // }
184 
185  double mult = (double)trackmap->size();
186  cout << " Number of tracks = " << trackmap->size() << endl;
187  unsigned int centbin = 99999;
188  for(unsigned int i=0; i<NCENT; i++) { if(mult>_multbins[i] && mult <=_multbins[i+1]) { centbin = i; } }
189  if(centbin<0 || centbin>=NCENT) {cout << "BAD MULTIPLICITY = " << mult << endl; return Fun4AllReturnCodes::ABORTEVENT; }
190  cout << "Centrality bin = " << centbin << endl;
191 
192  double _vtxbinsize = (_ZMAX - _ZMIN)/double(NZ);
193  //cout << " _vtxbinsize = " << _vtxbinsize << endl;
194 
195  vector<sPHElectronv1> elepos;
196 // vector<sPHElectronv1> electrons;
197 // vector<sPHElectronv1> positrons;
198 
199 // My own electron ID
200 /*
201  for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
202  {
203  SvtxTrack *track = iter->second;
204  if(!isElectron(track)) continue;
205  double px = track->get_px();
206  double py = track->get_py();
207  double pt = sqrt(px*px + py*py);
208  int charge = track->get_charge();
209  double x = track->get_x();
210  double y = track->get_y();
211  double z = track->get_z();
212  unsigned int vtxbin = (z - _ZMIN)/_vtxbinsize;
213  if(vtxbin<0 || vtxbin>=NZ) continue;
214  unsigned int vtxid = track->get_vertex_id();
215  if(vtxid<0 || vtxid>=global_vtxmap->size()) continue;
216  cout << "electron: "<<charge<<" "<<pt<<" "<<x<<" "<<y<<" "<<z<<" "<<vtxid<<" "<<vtxbin<< endl;
217  GlobalVertex* gvtx = global_vtxmap->get(vtxid);
218  cout << "global vertex: "<<gvtx->get_x()<<" "<<gvtx->get_y()<<" "<<gvtx->get_z()<<endl;
219  sPHElectronv1 tmpel = sPHElectronv1(track);
220  tmpel.set_zvtx(gvtx->get_z());
221  elepos.push_back(tmpel);
222  (_buffer[vtxbin][centbin]).push_back(tmpel);
223  } // end loop over tracks
224 */
225 
226  for(auto it = electrons.first; it != electrons.second; ++it)
227  {
228  SvtxTrack *track = trackmap->get(it->second);
229  unsigned int cemc_clusid = track->get_cal_cluster_id(SvtxTrack::CAL_LAYER::CEMC);
230  TrkrDefs::cluskey cemc_cluskey = track->get_cal_cluster_key(SvtxTrack::CAL_LAYER::CEMC);
231  float cemc_cluse = track->get_cal_cluster_e(SvtxTrack::CAL_LAYER::CEMC);
232  cout << "CEMC match: " << cemc_clusid << " " << cemc_cluskey << " " << cemc_cluse << endl;
233  RawCluster* cluster = cemc_clusters->getCluster(cemc_clusid);
234  double ee = cluster->get_energy();
235  double ecore = cluster->get_ecore();
236  double prob = cluster->get_prob();
237  double cemc_chi2 = cluster->get_chi2();
238  cout << "cluster: " << ee << " " << ecore << " " << prob << " " << cemc_chi2 << endl;
239 
240  double px = track->get_px();
241  double py = track->get_py();
242  double pt = sqrt(px*px + py*py);
243  int charge = track->get_charge();
244  double x = track->get_x();
245  double y = track->get_y();
246  double z = track->get_z();
247  unsigned int vtxbin = (z - _ZMIN)/_vtxbinsize;
248  if(vtxbin<0 || vtxbin>=NZ) continue;
249  unsigned int vtxid = track->get_vertex_id();
250  if(vtxid<0 || vtxid>=global_vtxmap->size()) continue;
251  cout << "electron: "<<charge<<" "<<pt<<" "<<x<<" "<<y<<" "<<z<<" "<<vtxid<<" "<<vtxbin<< endl;
252  GlobalVertex* gvtx = global_vtxmap->get(vtxid);
253  cout << "global vertex: "<<gvtx->get_x()<<" "<<gvtx->get_y()<<" "<<gvtx->get_z()<<endl;
254 
255  sPHElectronv1 tmpel = sPHElectronv1(track);
256  tmpel.set_zvtx(gvtx->get_z());
257  tmpel.set_cemc_ecore(ecore);
258  tmpel.set_cemc_prob(prob);
259  tmpel.set_cemc_chi2(cemc_chi2);
260 
261  elepos.push_back(tmpel);
262  (_buffer[vtxbin][centbin]).push_back(tmpel);
263 
264  }
265 
266  cout << "# of electrons/positrons = " << elepos.size() << endl;
267 
268  if(elepos.size()>1) {
269  for(long unsigned int i=0; i<elepos.size()-1; i++) {
270  for(long unsigned int j=i+1; j<elepos.size(); j++) {
271  sPHElectronPairv1 pair = sPHElectronPairv1(&elepos[i],&elepos[j]);
272  double mass = pair.get_mass();
273  int charge1 = (elepos[i]).get_charge();
274  int charge2 = (elepos[j]).get_charge();
275  int type = 0;
276  if(charge1*charge2<0) {type = 1;}
277  else if (charge1>0 && charge2>0) {type=2;}
278  else if (charge1<0 && charge2<0) {type=3;}
279  else {cout << "ERROR: wrong charge!" << endl;}
280  cout << "MASS = " << type << " " << mass << endl;
281  pair.set_type(type);
282  eePairs->insert(&pair);
283  }
284  }
285  }
286 
287 /*
288  for(long unsigned int i=0; i<electrons.size(); i++) {
289  for(long unsigned int j=0; j<positrons.size(); j++) {
290  sPHElectronPairv1 pair = sPHElectronPairv1(&electrons[i],&positrons[j]);
291  pair.set_type(1);
292  double mass = pair.get_mass();
293  cout << "MASS = " << mass << endl;
294  eePairs->insert(&pair);
295  }}
296 */
297 
298  int nmix = MakeMixedPairs(elepos, eePairs, centbin);
299  cout << "number of mixed pairs = " << nmix << endl;
300 
301  return Fun4AllReturnCodes::EVENT_OK;
302 }
303 
304 //======================================================================
305 
306 int PairMaker::MakeMixedPairs(std::vector<sPHElectronv1> elepos, sPHElectronPairContainerv1* eePairs, unsigned int centbin) {
307  int count=0;
308  double _vtxbinsize = (_ZMAX - _ZMIN)/double(NZ);
309 
310  for(unsigned int k=0; k<elepos.size(); k++) {
311 
312  sPHElectronv1 thisel = elepos[k];
313  double z = thisel.get_zvtx();
314  unsigned int vtxbin = (z - _ZMIN)/_vtxbinsize;
315  if(vtxbin<0 || vtxbin>=NZ) continue;
316 
317  unsigned int _num_mixes = 3;
318  unsigned int buffsize = (_buffer[vtxbin][centbin]).size();
319 
320  if(buffsize >= _min_buffer_depth) {
321  for(unsigned int i=0; i<_num_mixes; i++) {
322  double rnd = _rng->Uniform();
323  unsigned int irnd = rnd * buffsize;
324  sPHElectronv1 mixel = (_buffer[vtxbin][centbin]).at(irnd);
325 
326  sPHElectronPairv1 pair = sPHElectronPairv1(&thisel,&mixel);
327  int charge1 = thisel.get_charge();
328  int charge2 = mixel.get_charge();
329  int type = 0;
330  if(charge1*charge2<0) {type = 4;}
331  else if (charge1>0 && charge2>0) {type=5;}
332  else if (charge1<0 && charge2<0) {type=6;}
333  else {cout << "ERROR: wrong charge!" << endl;}
334  pair.set_type(type);
335  eePairs->insert(&pair);
336  cout << "Inserted MIXED pair with mass = " << type << " " << pair.get_mass() << endl;
337  count++;
338  } // end i loop
339  (_buffer[vtxbin][centbin]).push_back(elepos[k]); // keep filling buffer
340  }
341  else {
342  (_buffer[vtxbin][centbin]).push_back(elepos[k]); // keep filling buffer
343  }
344 
345  } //end k loop over elepos entries
346 
347  return count;
348 }
349 
350 //======================================================================
351 
352 bool PairMaker::isElectron(SvtxTrack* trk)
353 {
354  double px = trk->get_px();
355  double py = trk->get_py();
356  double pz = trk->get_pz();
357  double pt = sqrt(px*px+py*py);
358  double pp = sqrt(pt*pt+pz*pz);
359  double e3x3 = trk->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::CEMC);
360  if(pt<2.0) return false;
361  if(pp==0.) return false;
362  if(isnan(e3x3)) return false;
363  if(e3x3/pp<0.7) return false;
364  double chisq = trk->get_chisq();
365  double ndf = trk->get_ndf();
366  if((chisq/ndf)>10.) return false;
367  int nmvtx = 0;
368  int ntpc = 0;
369  for (SvtxTrack::ConstClusterKeyIter iter = trk->begin_cluster_keys();
370  iter != trk->end_cluster_keys();
371  ++iter)
372  {
373  TrkrDefs::cluskey cluster_key = *iter;
374  int trackerid = TrkrDefs::getTrkrId(cluster_key);
375  if(trackerid==0) nmvtx++;
376  if(trackerid==2) ntpc++;
377  }
378  if(nmvtx<1) return false;
379  if(ntpc<20) return false;
380 
381  return true;
382 }
383 
384 //==============================================================================
385 
386 int PairMaker::End(PHCompositeNode *topNode)
387 {
388  cout << "END: ====================================" << endl;
389  cout << "mixing buffers: " << endl;
390  for(unsigned int i=0; i<NZ; i++) {
391  for(unsigned int j=0; j<NCENT; j++) {
392  cout << i << " " << j << " " << (_buffer[i][j]).size() << endl;
393  (_buffer[i][j]).clear();
394  }
395  }
396  cout << "=========================================" << endl;
397 
398  return Fun4AllReturnCodes::EVENT_OK;
399 }
400 
401