EventMix/src/PairMaker.cc

0001
0002 #include "PairMaker.h"
0003
0004 #include <cstdlib>
0005 #include <cstdio>
0006 #include <iostream>
0007 #include <iomanip>
0008 #include <fstream>
0009
0010 #include "TFile.h"
0011 #include "TLorentzVector.h"
0012 #include "TRandom2.h"
0013
0014 #include <trackbase_historic/SvtxTrack.h>
0015 #include <trackbase_historic/SvtxTrackMap.h>
0016 #include <trackbase_historic/SvtxVertex.h>
0017 #include <trackbase_historic/SvtxVertexMap.h>
0018 #include <trackbase/TrkrDefs.h>
0019
0020 #include <globalvertex/GlobalVertexMap.h>
0021 #include <globalvertex/GlobalVertex.h>
0022
0023 #include <calobase/RawClusterContainer.h>
0024 #include <calobase/RawCluster.h>
0025 #include <calobase/RawClusterv1.h>
0026
0027 #include <phool/getClass.h>
0028 #include <phool/recoConsts.h>
0029 #include <phool/PHCompositeNode.h>
0030 #include <phool/PHIODataNode.h>
0031 #include <phool/PHNodeIterator.h>
0032 #include <phool/PHRandomSeed.h>
0033 #include <fun4all/Fun4AllReturnCodes.h>
0034
0035 #include "sPHElectron.h"
0036 #include "sPHElectronv1.h"
0037 #include "sPHElectronPair.h"
0038 #include "sPHElectronPairv1.h"
0039 #include "sPHElectronPairContainer.h"
0040 #include "sPHElectronPairContainerv1.h"
0041
0042 #include "trackpidassoc/TrackPidAssoc.h"
0043
0044 //#include <gsl/gsl_rng.h>
0045
0046 typedef PHIODataNode<PHObject> PHObjectNode_t;
0047
0048 using namespace std;
0049
0050 //==============================================================
0051
0052 PairMaker::PairMaker(const std::string &name, const std::string &filename) : SubsysReco(name)
0053 {
0054   _ZMIN = -15.;
0055   _ZMAX = 15.;
0056   _multbins.push_back(0.);
0057 //  _multbins.push_back(75.);
0058 //  _multbins.push_back(250.);
0059 //  _multbins.push_back(900.);
0060   _multbins.push_back(3000.);
0061   _multbins.push_back(9999.);
0062   _min_buffer_depth = 10;
0063   _max_buffer_depth = 50;
0064   outnodename = "ElectronPairs";
0065   _rng = nullptr;
0066   EventNumber=0;
0067 }
0068
0069 //==============================================================
0070
0071 int PairMaker::Init(PHCompositeNode *topNode)
0072 {
0073
0074   _rng = new TRandom2();
0075   _rng->SetSeed(0);
0076
0077   PHNodeIterator iter(topNode);
0078   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0079   if (!dstNode)
0080   {
0081     cerr << PHWHERE << " ERROR: Can not find DST node." << endl;
0082     return Fun4AllReturnCodes::ABORTEVENT;
0083   }
0084
0085   sPHElectronPairContainerv1* eePairs = new sPHElectronPairContainerv1();
0086
0087   PHCompositeNode *PairNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", outnodename));
0088   if (!PairNode)
0089   {
0090     PHObjectNode_t *PairNode = new PHIODataNode<PHObject>(eePairs,outnodename.c_str(),"PHObject");
0091     dstNode->addNode(PairNode);
0092     cout << PHWHERE << " INFO: added " << outnodename << endl;
0093   }
0094   else { cout << PHWHERE << " INFO: " << outnodename << " already exists." << endl; }
0095
0096   std::cout << "PairMaker::Init ended." << endl;
0097   return Fun4AllReturnCodes::EVENT_OK;
0098 }
0099
0100 //==============================================================
0101
0102 int PairMaker::InitRun(PHCompositeNode *topNode)
0103 {
0104   return Fun4AllReturnCodes::EVENT_OK;
0105 }
0106
0107 //==============================================================
0108
0109 int PairMaker::process_event(PHCompositeNode *topNode)
0110 {
0111   return process_event_test(topNode);
0112 }
0113
0114 //======================================================================
0115
0116 int PairMaker::process_event_test(PHCompositeNode *topNode) {
0117   EventNumber++;
0118
0119 //  int howoften = 1;
0120 //  if((EventNumber-1)%howoften==0) {
0121     cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
0122 //  }
0123
0124   if(EventNumber==1) topNode->print();
0125
0126   sPHElectronPairContainerv1 *eePairs = findNode::getClass<sPHElectronPairContainerv1>(topNode,outnodename.c_str());
0127   if(!eePairs) { cerr << outnodename << " not found!" << endl; }
0128     else { cout << "Found " << outnodename << " node." << endl; }
0129
0130   GlobalVertexMap *global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0131   if(!global_vtxmap) {
0132     cerr << PHWHERE << " ERROR: Can not find GlobalVertexMap node." << endl;
0133     return Fun4AllReturnCodes::ABORTEVENT;
0134   }
0135   //cout << "Number of GlobalVertexMap entries = " << global_vtxmap->size() << endl;
0136
0137   SvtxVertexMap *vtxmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0138   if(!vtxmap) {
0139       cout << "SvtxVertexMap node not found!" << endl;
0140       return Fun4AllReturnCodes::ABORTEVENT;
0141   }
0142   //cout << "Number of SvtxVertexMap entries = " << vtxmap->size() << endl;
0143
0144   SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0145   if(!trackmap) {
0146     cerr << PHWHERE << " ERROR: Can not find SvtxTrackMap node." << endl;
0147     return Fun4AllReturnCodes::ABORTEVENT;
0148   }
0149
0150   RawClusterContainer* cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0151   if(!cemc_clusters) {
0152     cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
0153     return Fun4AllReturnCodes::ABORTEVENT;
0154   }
0155   else { cout << "FOUND CLUSTER_CEMC node." << endl; }
0156
0157   int mycount = 0;
0158   RawClusterContainer::ConstRange cluster_range = cemc_clusters->getClusters();
0159   for (RawClusterContainer::ConstIterator cluster_iter = cluster_range.first; cluster_iter != cluster_range.second; cluster_iter++)
0160   {
0161     //RawCluster *cluster = cluster_iter->second;
0162     //double phi = cluster->get_phi();
0163     //double z = cluster->get_z();
0164     //double ee = cluster->get_energy();
0165     //int ntowers = cluster->getNTowers();
0166     //if(ee>2.) cout << "cluster: " << ee << " " << ntowers << " " << phi << " " << z << endl;
0167     mycount++;
0168   }
0169   cout << "Number of CEMC clusters = " << mycount << endl;
0170
0171   TrackPidAssoc *track_pid_assoc =  findNode::getClass<TrackPidAssoc>(topNode, "TrackPidAssoc");
0172   if(!track_pid_assoc) {
0173     cerr << PHWHERE << "ERROR: CAN NOT FIND TrackPidAssoc Node!" << endl;
0174     return Fun4AllReturnCodes::ABORTEVENT;
0175   }
0176   auto electrons = track_pid_assoc->getTracks(TrackPidAssoc::electron);
0177
0178 //    for(auto it = electrons.first; it != electrons.second; ++it)
0179 //    {
0180 //      SvtxTrack *tr = trackmap->get(it->second);
0181 //      double p = tr->get_p();
0182 //      std::cout << " pid " << it->first << " track ID " << it->second << " mom " << p << std::endl;
0183 //    }
0184
0185   double mult = (double)trackmap->size();
0186   cout << "   Number of tracks = " << trackmap->size() << endl;
0187   unsigned int centbin = 99999;
0188   for(unsigned int i=0; i<NCENT; i++) { if(mult>_multbins[i] && mult <=_multbins[i+1]) { centbin = i; } }
0189   if(centbin<0 || centbin>=NCENT) {cout << "BAD MULTIPLICITY = " << mult << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0190   cout << "Centrality bin = " << centbin << endl;
0191
0192   double _vtxbinsize  = (_ZMAX - _ZMIN)/double(NZ);
0193   //cout << " _vtxbinsize = " << _vtxbinsize << endl;
0194
0195   vector<sPHElectronv1> elepos;
0196 //  vector<sPHElectronv1> electrons;
0197 //  vector<sPHElectronv1> positrons;
0198
0199 // My own electron ID
0200 /*
0201   for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
0202   {
0203     SvtxTrack *track = iter->second;
0204     if(!isElectron(track)) continue;
0205     double px = track->get_px();
0206     double py = track->get_py();
0207     double pt = sqrt(px*px + py*py);
0208     int charge = track->get_charge();
0209     double x = track->get_x();
0210     double y = track->get_y();
0211     double z = track->get_z();
0212     unsigned int vtxbin = (z - _ZMIN)/_vtxbinsize;
0213     if(vtxbin<0 || vtxbin>=NZ) continue;
0214     unsigned int vtxid = track->get_vertex_id();
0215     if(vtxid<0 || vtxid>=global_vtxmap->size()) continue;
0216     cout << "electron: "<<charge<<" "<<pt<<" "<<x<<" "<<y<<" "<<z<<" "<<vtxid<<" "<<vtxbin<< endl;
0217     GlobalVertex* gvtx = global_vtxmap->get(vtxid);
0218     cout << "global vertex: "<<gvtx->get_x()<<" "<<gvtx->get_y()<<" "<<gvtx->get_z()<<endl;
0219     sPHElectronv1 tmpel = sPHElectronv1(track);
0220     tmpel.set_zvtx(gvtx->get_z());
0221     elepos.push_back(tmpel);
0222     (_buffer[vtxbin][centbin]).push_back(tmpel);
0223   } // end loop over tracks
0224 */
0225
0226     for(auto it = electrons.first; it != electrons.second; ++it)
0227     {
0228       SvtxTrack *track = trackmap->get(it->second);
0229       unsigned int cemc_clusid = track->get_cal_cluster_id(SvtxTrack::CAL_LAYER::CEMC);
0230       TrkrDefs::cluskey cemc_cluskey = track->get_cal_cluster_key(SvtxTrack::CAL_LAYER::CEMC);
0231       float cemc_cluse = track->get_cal_cluster_e(SvtxTrack::CAL_LAYER::CEMC);
0232         cout << "CEMC match: " << cemc_clusid << " " << cemc_cluskey << " " << cemc_cluse << endl;
0233         RawCluster* cluster = cemc_clusters->getCluster(cemc_clusid);
0234         double ee = cluster->get_energy();
0235         double ecore = cluster->get_ecore();
0236         double prob = cluster->get_prob();
0237         double cemc_chi2 = cluster->get_chi2();
0238         cout << "cluster: " << ee << " " << ecore << " " << prob << " " << cemc_chi2 << endl;
0239
0240       double px = track->get_px();
0241       double py = track->get_py();
0242       double pt = sqrt(px*px + py*py);
0243       int charge = track->get_charge();
0244       double x = track->get_x();
0245       double y = track->get_y();
0246       double z = track->get_z();
0247       unsigned int vtxbin = (z - _ZMIN)/_vtxbinsize;
0248       if(vtxbin<0 || vtxbin>=NZ) continue;
0249       unsigned int vtxid = track->get_vertex_id();
0250       if(vtxid<0 || vtxid>=global_vtxmap->size()) continue;
0251       cout << "electron: "<<charge<<" "<<pt<<" "<<x<<" "<<y<<" "<<z<<" "<<vtxid<<" "<<vtxbin<< endl;
0252       GlobalVertex* gvtx = global_vtxmap->get(vtxid);
0253       cout << "global vertex: "<<gvtx->get_x()<<" "<<gvtx->get_y()<<" "<<gvtx->get_z()<<endl;
0254
0255       sPHElectronv1 tmpel = sPHElectronv1(track);
0256       tmpel.set_zvtx(gvtx->get_z());
0257       tmpel.set_cemc_ecore(ecore);
0258       tmpel.set_cemc_prob(prob);
0259       tmpel.set_cemc_chi2(cemc_chi2);
0260
0261       elepos.push_back(tmpel);
0262       (_buffer[vtxbin][centbin]).push_back(tmpel);
0263
0264     }
0265
0266   cout << "# of electrons/positrons = " << elepos.size() << endl;
0267
0268   if(elepos.size()>1) {
0269     for(long unsigned int i=0; i<elepos.size()-1; i++) {
0270       for(long unsigned int j=i+1; j<elepos.size(); j++) {
0271         sPHElectronPairv1 pair = sPHElectronPairv1(&elepos[i],&elepos[j]);
0272         double mass = pair.get_mass();
0273         int charge1 = (elepos[i]).get_charge();
0274         int charge2 = (elepos[j]).get_charge();
0275         int type = 0;
0276         if(charge1*charge2<0) {type = 1;}
0277           else if (charge1>0 && charge2>0) {type=2;}
0278             else if (charge1<0 && charge2<0) {type=3;}
0279               else {cout << "ERROR: wrong charge!" << endl;}
0280         cout << "MASS = " << type << " " << mass << endl;
0281         pair.set_type(type);
0282         eePairs->insert(&pair);
0283       }
0284     }
0285   }
0286
0287 /*
0288   for(long unsigned int i=0; i<electrons.size(); i++) {
0289   for(long unsigned int j=0; j<positrons.size(); j++) {
0290     sPHElectronPairv1 pair = sPHElectronPairv1(&electrons[i],&positrons[j]);
0291     pair.set_type(1);
0292     double mass = pair.get_mass();
0293     cout << "MASS = " << mass << endl;
0294     eePairs->insert(&pair);
0295   }}
0296 */
0297
0298   int nmix = MakeMixedPairs(elepos, eePairs, centbin);
0299   cout << "number of mixed pairs = " << nmix << endl;
0300
0301   return Fun4AllReturnCodes::EVENT_OK;
0302 }
0303
0304 //======================================================================
0305
0306 int PairMaker::MakeMixedPairs(std::vector<sPHElectronv1> elepos, sPHElectronPairContainerv1* eePairs, unsigned int centbin) {
0307   int count=0;
0308   double _vtxbinsize  = (_ZMAX - _ZMIN)/double(NZ);
0309
0310   for(unsigned int k=0; k<elepos.size(); k++) {
0311
0312     sPHElectronv1 thisel = elepos[k];
0313     double z = thisel.get_zvtx();
0314     unsigned int vtxbin = (z - _ZMIN)/_vtxbinsize;
0315     if(vtxbin<0 || vtxbin>=NZ) continue;
0316
0317     unsigned int _num_mixes = 3;
0318     unsigned int buffsize = (_buffer[vtxbin][centbin]).size();
0319
0320     if(buffsize >= _min_buffer_depth) {
0321       for(unsigned int i=0; i<_num_mixes; i++) {
0322         double rnd = _rng->Uniform();
0323         unsigned int irnd = rnd * buffsize;
0324         sPHElectronv1 mixel = (_buffer[vtxbin][centbin]).at(irnd);
0325
0326         sPHElectronPairv1 pair = sPHElectronPairv1(&thisel,&mixel);
0327         int charge1 = thisel.get_charge();
0328         int charge2 = mixel.get_charge();
0329         int type = 0;
0330         if(charge1*charge2<0) {type = 4;}
0331           else if (charge1>0 && charge2>0) {type=5;}
0332             else if (charge1<0 && charge2<0) {type=6;}
0333               else {cout << "ERROR: wrong charge!" << endl;}
0334         pair.set_type(type);
0335         eePairs->insert(&pair);
0336         cout << "Inserted MIXED pair with mass = " << type << " " << pair.get_mass() << endl;
0337         count++;
0338       } // end i loop
0339       (_buffer[vtxbin][centbin]).push_back(elepos[k]);  // keep filling buffer
0340     }
0341     else {
0342       (_buffer[vtxbin][centbin]).push_back(elepos[k]);  // keep filling buffer
0343     }
0344
0345   } //end k loop over elepos entries
0346
0347   return count;
0348 }
0349
0350 //======================================================================
0351
0352 bool PairMaker::isElectron(SvtxTrack* trk)
0353 {
0354   double px = trk->get_px();
0355   double py = trk->get_py();
0356   double pz = trk->get_pz();
0357   double pt = sqrt(px*px+py*py);
0358   double pp = sqrt(pt*pt+pz*pz);
0359   double e3x3 = trk->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::CEMC);
0360   if(pt<2.0) return false;
0361   if(pp==0.) return false;
0362   if(isnan(e3x3)) return false;
0363   if(e3x3/pp<0.7) return false;
0364   double chisq = trk->get_chisq();
0365   double ndf = trk->get_ndf();
0366   if((chisq/ndf)>10.) return false;
0367   int nmvtx = 0;
0368   int ntpc = 0;
0369   for (SvtxTrack::ConstClusterKeyIter iter = trk->begin_cluster_keys();
0370      iter != trk->end_cluster_keys();
0371      ++iter)
0372   {
0373     TrkrDefs::cluskey cluster_key = *iter;
0374     int trackerid = TrkrDefs::getTrkrId(cluster_key);
0375     if(trackerid==0) nmvtx++;
0376     if(trackerid==2) ntpc++;
0377   }
0378   if(nmvtx<1) return false;
0379   if(ntpc<20) return false;
0380
0381   return true;
0382 }
0383
0384 //==============================================================================
0385
0386 int PairMaker::End(PHCompositeNode *topNode)
0387 {
0388   cout << "END: ====================================" << endl;
0389   cout << "mixing buffers: " << endl;
0390   for(unsigned int i=0; i<NZ; i++) {
0391     for(unsigned int j=0; j<NCENT; j++) {
0392       cout << i << " " << j << "    " << (_buffer[i][j]).size() << endl;
0393       (_buffer[i][j]).clear();
0394     }
0395   }
0396   cout << "=========================================" << endl;
0397
0398   return Fun4AllReturnCodes::EVENT_OK;
0399 }
0400
0401