Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:12:34 

0001 
0002 #include "sPHAnalysis.h"
0003 
0004 #include <cmath>
0005 #include <cstdlib>
0006 #include <cstdio>
0007 #include <iostream>
0008 #include <iomanip>
0009 #include <fstream>
0010 
0011 #include "TFile.h"
0012 #include "TNtuple.h"
0013 #include "TH1D.h"
0014 #include "TF1.h"
0015 #include "TLorentzVector.h"
0016 #include "TRandom2.h"
0017 
0018 #include <trackbase_historic/SvtxTrack.h>
0019 #include <trackbase_historic/SvtxTrackMap.h>
0020 #include <trackbase_historic/SvtxVertex.h>
0021 #include <trackbase_historic/SvtxVertexMap.h>
0022 #include <trackbase_historic/TrackSeed.h>
0023 #include <trackbase_historic/TrackSeed_v1.h>
0024 #include <trackbase_historic/SvtxTrackSeed_v1.h>
0025 #include <trackbase_historic/TrackSeedContainer.h>
0026 #include <trackbase/TrkrCluster.h>
0027 #include <trackbase/TrkrClusterv4.h>
0028 #include <trackbase/TrkrClusterContainer.h>
0029 #include <trackbase/TrkrDefs.h>
0030 
0031 #include <g4vertex/GlobalVertexMap.h>
0032 #include <g4vertex/GlobalVertex.h>
0033 
0034 #include <calobase/RawClusterContainer.h>
0035 #include <calobase/RawCluster.h>
0036 #include <calobase/RawClusterv1.h>
0037 #include <calobase/RawTowerv2.h>
0038 #include <calobase/RawTowerGeom.h>
0039 #include <calobase/RawTowerContainer.h>
0040 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0041 #include <calobase/RawTowerGeomContainer.h>
0042 
0043 #include <phhepmc/PHHepMCGenEvent.h>
0044 #include <phhepmc/PHHepMCGenEventMap.h>
0045 
0046 //#include <HepMC/GenEvent.h>              // for GenEvent::particle_const_ite...
0047 //#include <HepMC/GenParticle.h>           // for GenParticle
0048 #include <HepMC/GenVertex.h>             // for GenVertex, GenVertex::partic...
0049 #include <HepMC/IteratorRange.h>         // for ancestors, children, descend...
0050 #include <HepMC/SimpleVector.h>   // for FourVector
0051 
0052 #include <phool/getClass.h>
0053 #include <phool/recoConsts.h>
0054 #include <phool/PHCompositeNode.h>
0055 #include <phool/PHIODataNode.h>
0056 #include <phool/PHRandomSeed.h>
0057 #include <fun4all/Fun4AllReturnCodes.h>
0058 
0059 #include <g4main/PHG4TruthInfoContainer.h>
0060 #include <g4main/PHG4Particle.h>
0061 #include <g4main/PHG4VtxPoint.h>
0062 
0063 #include <gsl/gsl_rng.h>
0064 
0065 #include "/gpfs/mnt/gpfs02/sphenix/user/lebedev/mdc/analysis/install/include/eventmix/sPHElectron.h"
0066 #include "/gpfs/mnt/gpfs02/sphenix/user/lebedev/mdc/analysis/install/include/eventmix/sPHElectronv1.h"
0067 #include "/gpfs/mnt/gpfs02/sphenix/user/lebedev/mdc/analysis/install/include/eventmix/sPHElectronPair.h"
0068 #include "/gpfs/mnt/gpfs02/sphenix/user/lebedev/mdc/analysis/install/include/eventmix/sPHElectronPairv1.h"
0069 #include "/gpfs/mnt/gpfs02/sphenix/user/lebedev/mdc/analysis/install/include/eventmix/sPHElectronPairContainer.h"
0070 #include "/gpfs/mnt/gpfs02/sphenix/user/lebedev/mdc/analysis/install/include/eventmix/sPHElectronPairContainerv1.h"
0071 
0072 using namespace std;
0073 //using namespace HepMC;
0074 
0075 //ofstream ofs;
0076 //std::ofstream ascii_io( "text_hepmc.txt" );
0077 
0078 //==============================================================
0079 
0080 sPHAnalysis::sPHAnalysis(const std::string &name, const std::string &filename) : SubsysReco(name)
0081 {
0082   OutputNtupleFile=NULL;
0083   OutputFileName=filename;
0084   EventNumber=0;
0085   ntp1=NULL;
0086   hdeta=NULL;
0087   hdphi=NULL;
0088   hmass=NULL;
0089   hgmass=NULL;
0090   hgmass0=NULL;
0091   hgmass09=NULL;
0092   heop=NULL;
0093   heop3x3=NULL;
0094   heop5x5=NULL;
0095   hdca2d=NULL;
0096   hdca3dxy=NULL;
0097   hdca3dz=NULL;
0098   hchi2=NULL;
0099   hndf=NULL;
0100   hquality=NULL;
0101 
0102   _whattodo = 2;
0103   _rng = nullptr;
0104 }
0105 
0106 //==============================================================
0107 
0108 int sPHAnalysis::Init(PHCompositeNode *topNode) 
0109 {
0110 
0111 //  ofs.open("test.txt");
0112 
0113   std::cout << "sPHAnalysis::Init started..." << endl;
0114   OutputNtupleFile = new TFile(OutputFileName.c_str(),"RECREATE");
0115   std::cout << "sPHAnalysis::Init: output file " << OutputFileName.c_str() << " opened." << endl;
0116 
0117   ntppid = new  TNtuple("ntppid","","charge:pt:eta:mom:nmvtx:ntpc:chi2:ndf:cemc_ecore:cemc_e3x3:cemc_prob:cemc_chi2:cemc_dphi:cemc_deta:quality:dca2d:hcalin_e3x3:hcalout_e3x3:gdist:cemc-dphi_3x3:cemc_deta_3x3:hcalin_dphi:hcalin_deta:hcalout_dphi:hcalout_deta:hcalin_clusdphi:hcalin_clusdeta:hcalin_cluse:hcalout_clusdphi:hcalout_clusdeta:hcalout_cluse");
0118 
0119   ntp1 = new  TNtuple("ntp1","","type:mass:pt:eta:pt1:pt2:e3x31:e3x32:emce1:emce2:p1:p2:chisq1:chisq2:dca2d1:dca2d2:dca3dxy1:dca3dxy2:dca3dz1:dcz3dz2:mult:rap:nmvtx1:nmvtx2:ntpc1:ntpc2");
0120 
0121   ntpmc1 = new TNtuple("ntpmc1","","charge:pt:eta");
0122   ntpmc2 = new TNtuple("ntpmc2","","type:mass:pt:eta:rap:pt1:pt2:eta1:eta2");
0123 
0124   ntp2   = new TNtuple("ntp2","",  "type:mass:pt:eta:rap:pt1:pt2:eta1:eta2:mult:emult:pmult:chisq1:chisq2:dca2d1:dca2d2:eop3x3_1:eop3x3_2:dretaphi:nmvtx1:nmvtx2:nintt1:nintt2:ntpc1:ntpc2:b");
0125 
0126   ntpb = new TNtuple("ntpb","","bimp:mult:truemult:cemcmult:evtno:evtno5:ginacc:ngood");
0127 
0128   ntp_notracking = new TNtuple("ntp_notracking","","mass:pt:pt1:pt2:hoe1:hoe2");
0129   h_notracking_etabins_em =  new TH1D("h_notracking_etabins_em","",200,-0.5,199.5);
0130   h_notracking_phibins_em =  new TH1D("h_notracking_phibins_em","",400,-0.5,399.5);
0131   h_notracking_etabins =  new TH1D("h_notracking_etabins","",100,-0.5,99.5);
0132   h_notracking_phibins =  new TH1D("h_notracking_phibins","",100,-0.5,99.5);
0133   h_notracking_etabinsout =  new TH1D("h_notracking_etabinsout","",100,-0.5,99.5);
0134   h_notracking_phibinsout =  new TH1D("h_notracking_phibinsout","",100,-0.5,99.5);
0135   h_notracking_eoh =  new TH1D("h_notracking_eoh","",10000,0.0,1.0);
0136 
0137   hmult =  new TH1D("hmult","",100,0.,2000.);
0138   hmass =  new TH1D("hmass","",160,4.,12.);
0139   hgmass =  new TH1D("hgmass","",160,4.,12.);
0140   hgmass0 =  new TH1D("hgmass0","",160,4.,12.);
0141   hgmass09 =  new TH1D("hgmass09","",160,4.,12.);
0142   heop3x3 =  new TH1D("heop3x3","",180,0.,1.8);
0143   heop5x5 =  new TH1D("heop5x5","",180,0.,1.8);
0144   heop =  new TH1D("heop","",180,0.,1.8);
0145   hdphi =  new TH1D("hdphi","",200,0.25,0.25);
0146   hdeta = new TH1D("hdeta","",200,0.25,0.25);
0147 
0148   hdca2d = new TH1D("hdca2d","",1000,-0.1,0.1);
0149   hdca3dxy = new TH1D("hdca3dxy","",1000,-0.1,0.1);
0150   hdca3dz = new TH1D("hdca3dz","",1000,-0.1,0.1);
0151   hchi2 = new TH1D("hchi2","",1000,0.,50.);
0152   hndf = new TH1D("hndf","",1000,0.,50.);
0153   hquality = new TH1D("hquality","",1000,0.,50.);
0154 
0155   test = new TH1D("test","",1000,0.,50.);
0156 
0157 //  HepMC::write_HepMC_IO_block_begin( ascii_io );
0158 
0159   _rng = new TRandom2();
0160   _rng->SetSeed(0);
0161 
0162 //  fsin = new TF1("fsin", "sin(x)", 0, M_PI);
0163 
0164   std::cout << "sPHAnalysis::Init ended." << endl;
0165   return Fun4AllReturnCodes::EVENT_OK;
0166 }
0167 
0168 //==============================================================
0169   
0170 int sPHAnalysis::InitRun(PHCompositeNode *topNode) 
0171 {
0172   std::cout << "sPHAnalysis::InitRun started..." << std::endl;
0173   std::cout << "sPHAnalysis::InitRun ended." << std::endl;
0174   return Fun4AllReturnCodes::EVENT_OK;
0175 }
0176 
0177 //==============================================================
0178 
0179 int sPHAnalysis::process_event(PHCompositeNode *topNode) 
0180 {
0181   //return process_event_bimp(topNode);
0182   //return process_event_hepmc(topNode);
0183   if(_whattodo==0) {
0184 //    cout << "MAKING PAIRS." << endl;
0185     return process_event_pairs(topNode);
0186   } else if(_whattodo==1) {
0187 //    cout << "DOING ANALYSIS." << endl;
0188     return process_event_test(topNode);
0189   } else if(_whattodo==2) {
0190 //    cout << "PYTHIA UPSILONS." << endl;
0191     return process_event_pythiaupsilon(topNode);
0192   } else if(_whattodo==3) {
0193         return process_event_upsilons(topNode);
0194   } else if(_whattodo==4) {
0195       return process_event_notracking(topNode);
0196   } else if(_whattodo==5) {
0197       return process_event_filtered(topNode);
0198   } else { cerr << "ERROR: wrong choice of what to do." << endl; return Fun4AllReturnCodes::ABORTRUN; }
0199 }
0200 
0201 //======================================================================
0202 
0203 int sPHAnalysis::process_event_hepmc(PHCompositeNode *topNode) {
0204   EventNumber++;
0205 
0206   cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
0207   if(EventNumber==1) topNode->print();
0208 
0209   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0210     PHHepMCGenEvent *genevt = nullptr;
0211 
0212   for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
0213   {
0214     genevt = iter->second;
0215     if(!genevt) {cout<<"no phgenevt!" << endl; return Fun4AllReturnCodes::ABORTEVENT;}
0216     cout << "got PHHepMCGenEvent... " << genevt << endl;
0217   }
0218 
0219     HepMC::GenEvent *event = genevt->getEvent();
0220     if (!event) { cout << PHWHERE << "no genevent!" << endl; return Fun4AllReturnCodes::ABORTEVENT;}
0221     int npart = event->particles_size();
0222     int nvert = event->vertices_size();
0223     cout << "Number of particles = " << npart << " " << nvert << endl;
0224 
0225 /*
0226     for (HepMC::GenEvent::particle_const_iterator p = event->particles_begin(); p != event->particles_end(); ++p)
0227     {
0228       int pid = (*p)->pdg_id();
0229       int status = (*p)->status();
0230       double pt = ((*p)->momentum()).perp();
0231 //      double mass  = ((*p)->momentum()).m();
0232 //      double eta  = ((*p)->momentum()).eta();
0233       if(pt>2.0) cout << pid << " " << pt << " " << status << endl;
0234     }
0235 */
0236 
0237 // generate Upsilon
0238 
0239    double pt = _rng->Uniform() * 10.;
0240    double y  = (_rng->Uniform()*2.) - 1.;
0241    double phi = (2.0 * M_PI) * _rng->Uniform();
0242    double mass = 9.3987;
0243    double mt = sqrt(mass * mass + pt * pt);
0244    double eta = asinh(sinh(y) * mt / pt);
0245    TLorentzVector vm;
0246    vm.SetPtEtaPhiM(pt, eta, phi, mass);
0247 
0248 // decay it to electrons
0249 
0250     double m1 = 0.000511;
0251     double m2 = 0.000511;
0252 
0253 // first electroin in Upsilon rest frame
0254     double E1 = (mass * mass - m2 * m2 + m1 * m1) / (2.0 * mass);
0255     double p1 = sqrt((mass * mass - (m1 + m2) * (m1 + m2)) * (mass * mass - (m1 - m2) * (m1 - m2))) / (2.0 * mass);
0256     double phi1 = _rng->Uniform() * 2. * M_PI;
0257     double th1  = fsin->GetRandom();
0258     double px1 = p1 * sin(th1) * cos(phi1);
0259     double py1 = p1 * sin(th1) * sin(phi1);
0260     double pz1 = p1 * cos(th1);
0261     TLorentzVector v1;
0262     v1.SetPxPyPzE(px1, py1, pz1, E1);
0263 
0264 // boost to lab frame
0265     double betax = vm.Px() / vm.E();
0266     double betay = vm.Py() / vm.E();
0267     double betaz = vm.Pz() / vm.E();
0268     v1.Boost(betax, betay, betaz);
0269 
0270 // second electron
0271     TLorentzVector v2 = vm - v1;    
0272     cout << "decay electrons: " << v1.Pt() << " " << v2.Pt() << endl;
0273 
0274     HepMC::FourVector fv1 = HepMC::FourVector(v1.Px(),v1.Py(),v1.Pz());
0275     HepMC::GenParticle gp1 = HepMC::GenParticle(fv1,-11);
0276     HepMC::FourVector fv2 = HepMC::FourVector(v2.Px(),v2.Py(),v2.Pz());
0277     HepMC::GenParticle gp2 = HepMC::GenParticle(fv2,11);
0278 
0279     HepMC::GenEvent* newevent = new HepMC::GenEvent(*event);
0280 
0281     for (HepMC::GenEvent::vertex_iterator v = newevent->vertices_begin(); v != newevent->vertices_end(); ++v)
0282     {
0283        int npin = (*v)->particles_in_size();
0284        int npout = (*v)->particles_out_size();
0285        HepMC::ThreeVector pnt = (*v)->point3d();
0286        cout << "event vertex: " << npin << " " << npout << " " << pnt.x() << " " << pnt.y() << " " << pnt.z() << endl;
0287        cout << "adding particle with id " << gp1.pdg_id() << endl;
0288        (*v)->add_particle_out(&gp1);
0289        cout << "   after adding electron1: " << (*v)->particles_out_size() << endl;
0290        (*v)->add_particle_out(&gp2);
0291        cout << "   after adding electron2: " << (*v)->particles_out_size() << endl;
0292        break;  // add electrons to the first vertex
0293     }
0294 
0295     PHHepMCGenEvent* newgenevt = new PHHepMCGenEvent();
0296     newgenevt->addEvent(newevent);
0297 
0298 //    genevt->clearEvent();
0299     PHHepMCGenEvent* tmpptr = genevtmap->insert_event(1, newgenevt);
0300     cout << "inserted event " << tmpptr << endl;
0301 
0302   cout << "procedd_event_hepmc() ended." << endl;
0303   return Fun4AllReturnCodes::EVENT_OK;
0304 }
0305 
0306 //======================================================================
0307 
0308 int sPHAnalysis::process_event_pairs(PHCompositeNode *topNode) {
0309   EventNumber++;
0310   float tmp[99];
0311   cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
0312     //if(EventNumber==1) topNode->print();
0313 
0314   sPHElectronPairContainerv1 *eePairs = findNode::getClass<sPHElectronPairContainerv1>(topNode,"ElectronPairs");
0315   cout << eePairs << endl;
0316   if(!eePairs) { cout << "ERROR: ElectronPairs node not found!" << endl; }
0317     else { cout << "Found ElectronPairs node." << endl; }
0318 
0319   float mult = eePairs->size();
0320   cout << "Number of pairs = " << eePairs->size() << endl;
0321 
0322   for (sPHElectronPairContainerv1::ConstIter iter = eePairs->begin(); iter != eePairs->end(); ++iter)
0323   {
0324     sPHElectronPairv1* pair = iter->second;
0325     sPHElectronv1* electron = (sPHElectronv1*)pair->get_first();
0326     sPHElectronv1* positron = (sPHElectronv1*)pair->get_second();
0327     tmp[0] = pair->get_type();
0328     tmp[1] = pair->get_mass();
0329     tmp[2] = pair->get_pt();
0330     tmp[3] = pair->get_eta();
0331       double px1 = electron->get_px();
0332       double py1 = electron->get_py();
0333       double pz1 = electron->get_pz();
0334       double px2 = positron->get_px();
0335       double py2 = positron->get_py();
0336       double pz2 = positron->get_pz();
0337     tmp[4] = sqrt(px1*px1+py1*py1);
0338     tmp[5] = sqrt(px2*px2+py2*py2);
0339     tmp[6] = electron->get_e3x3();
0340     tmp[7] = positron->get_e3x3();
0341     tmp[8] = electron->get_emce();
0342     tmp[9] = positron->get_emce();
0343     tmp[10] = sqrt(px1*px1+py1*py1+pz1*pz1);
0344     tmp[11] = sqrt(px2*px2+py2*py2+pz2*pz2);
0345       double chisq1 = electron->get_chi2();
0346       double ndf1 = (double)electron->get_ndf();
0347       if(ndf1!=0.) { chisq1 = chisq1/ndf1;} else {chisq1 = 99999.;}
0348       double chisq2 = positron->get_chi2();
0349       double ndf2 = (double)positron->get_ndf();
0350       if(ndf2!=0.) { chisq2 = chisq2/ndf2;} else {chisq2 = 99999.;}
0351     tmp[12] = chisq1;
0352     tmp[13] = chisq2;
0353     tmp[14] = electron->get_dca2d();
0354     tmp[15] = positron->get_dca2d();
0355     tmp[16] = electron->get_dca3d_xy();
0356     tmp[17] = positron->get_dca3d_xy();
0357     tmp[18] = electron->get_dca3d_z();
0358     tmp[19] = positron->get_dca3d_z();
0359     tmp[20] = mult;
0360     tmp[21] = 0.; // rapidity
0361     tmp[22] = positron->get_nmvtx();
0362     tmp[23] = electron->get_nmvtx();
0363     tmp[24] = positron->get_ntpc();
0364     tmp[25] = electron->get_ntpc();
0365 
0366        ntp1->Fill(tmp);
0367   }
0368 
0369 
0370   cout << "process_event_pairs() ended." << endl;
0371   return Fun4AllReturnCodes::EVENT_OK;
0372 }
0373 
0374 //============================================================================
0375 
0376 int sPHAnalysis::process_event_bimp(PHCompositeNode *topNode) {
0377   EventNumber++;
0378   if((EventNumber-1)%10==0) {cout<<"------------------ EventNumber = " << EventNumber-1 << endl;}
0379   float tmp[99];
0380 
0381   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0382   PHHepMCGenEvent *genevt = nullptr;
0383 
0384   for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
0385   {
0386     genevt = iter->second;
0387     if(!genevt) {cout<<"ERROR: no PHHepMCGenEvent!" << endl; return Fun4AllReturnCodes::ABORTEVENT;}
0388   }
0389 
0390     HepMC::GenEvent *event = genevt->getEvent();
0391     if (!event) { cout << PHWHERE << "ERROR: no HepMC::GenEvent!" << endl; return Fun4AllReturnCodes::ABORTEVENT;}
0392 
0393     HepMC::HeavyIon* hi = event->heavy_ion();
0394     float impact_parameter = hi->impact_parameter();
0395     //cout << "HepMC::GenEvent: impact parameter = " << impact_parameter << endl;
0396 
0397     tmp[0] = impact_parameter;
0398     ntpb->Fill(tmp);
0399 
0400   return Fun4AllReturnCodes::EVENT_OK;
0401 }
0402 
0403 //=============================================================================================
0404 
0405 HepMC::GenParticle*  sPHAnalysis::GetParent(HepMC::GenParticle* p, HepMC::GenEvent* event)
0406 {
0407 
0408   HepMC::GenParticle* parent = NULL;
0409 //  if(!p->production_vertex()) return parent;
0410 
0411   for ( HepMC::GenVertex::particle_iterator mother = p->production_vertex()-> particles_begin(HepMC::ancestors);
0412         mother != p->production_vertex()-> particles_end(HepMC::ancestors);
0413         ++mother )
0414     {
0415           //if(abs((*mother)->pdg_id()) == 23 || abs((*mother)->pdg_id()) == 21) // Z0 or gluon
0416           //cout << "      mother pid = " << (*mother)->pdg_id() << endl;
0417           if(abs((*mother)->pdg_id()) == 553) // Upsilon
0418             {
0419               parent = *mother;
0420               break;
0421             }
0422       if(parent != NULL) break;
0423     }
0424 
0425   return parent;
0426 }
0427 
0428 
0429 //======================================================================
0430 
0431 int sPHAnalysis::process_event_pythiaupsilon(PHCompositeNode *topNode) {
0432   EventNumber++;
0433   int howoften = 1; 
0434   if((EventNumber-1)%howoften==0) { cout<<"------------ EventNumber = " << EventNumber-1 << endl; }
0435   if(EventNumber==1) topNode->print();
0436 
0437   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0438   PHHepMCGenEvent *genevt = nullptr;
0439 
0440   for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
0441   {
0442     genevt = iter->second;
0443     if(!genevt) {cout<<"ERROR: no PHHepMCGenEvent!" << endl; return Fun4AllReturnCodes::ABORTEVENT;}
0444     //cout << "got PHHepMCGenEvent... " << genevt << endl;
0445   }
0446     
0447     HepMC::GenEvent *event = genevt->getEvent();
0448     if (!event) { cout << PHWHERE << "ERROR: no HepMC::GenEvent!" << endl; return Fun4AllReturnCodes::ABORTEVENT;}
0449     
0450     int npart = event->particles_size();
0451     int nvert = event->vertices_size();
0452     cout << "HepMC::GenEvent: Number of particles, vertuces = " << npart << " " << nvert << endl;
0453     test->Fill(npart);
0454     
0455     bool einacc = false;
0456     bool pinacc = false;
0457     bool einacc2 = false;
0458     bool pinacc2 = false;
0459     double jpsipt = 0.;
0460     double jpsieta = 0.;
0461     //double jpsimass = 0.;
0462     double jpsirap = 0.;
0463     double pt1 = 0.;
0464     double pt2 = 0.;
0465     double eta1 = 0.;
0466     double eta2 = 0.;
0467     HepMC::GenVertex* jpsidecay_vtx = NULL;
0468 
0469 // Find Upsilon or J/psi
0470     for (HepMC::GenEvent::particle_const_iterator p = event->particles_begin(); p != event->particles_end(); ++p)
0471     {
0472       int pid = (*p)->pdg_id();
0473       int status = (*p)->status();
0474       double pt = ((*p)->momentum()).perp();
0475       double pz = ((*p)->momentum()).pz();
0476       double mass  = ((*p)->momentum()).m();
0477       double eta  = ((*p)->momentum()).eta();
0478       double ee = ((*p)->momentum()).e();
0479       double rap = 0.5*TMath::Log((ee+pz)/(ee-pz));
0480       
0481       if(abs(pid)==553 && status==2) {
0482         cout << "Upsilon: " << pid << " " << status << " " << pt << " " << mass << " " << eta << endl;
0483         HepMC::GenVertex* upsvtx = (*p)->production_vertex();
0484         cout << "   Upsilon production Z vertex = " << upsvtx->point3d().z() << endl; 
0485       }
0486 
0487       if(abs(pid)==443 && status==2) {
0488         cout << "J/psi: " << pid << " " << status << " " << pt << " " << mass << " " << eta << endl;
0489         jpsidecay_vtx = (*p)->end_vertex();
0490         jpsipt = pt;
0491         jpsieta = eta;
0492         jpsirap = rap;
0493       }
0494 
0495     }
0496 
0497 // Find decay electrons
0498     for (HepMC::GenEvent::particle_const_iterator p = event->particles_begin(); p != event->particles_end(); ++p)
0499     {
0500       int pid = (*p)->pdg_id();
0501       int status = (*p)->status();
0502       double pt = ((*p)->momentum()).perp();
0503 //      double pz = ((*p)->momentum()).pz();
0504       double mass  = ((*p)->momentum()).m();
0505       double eta  = ((*p)->momentum()).eta();
0506 //      double ee = ((*p)->momentum()).e();
0507 //      double rap = 0.5*ln((ee+pz)/(ee-pz));
0508 
0509       if(pid==-11  && status==1 && pt>0.0) {
0510         //HepMC::GenParticle* parent = GetParent(*p, event);
0511         //int parentid = 0; if(parent) {parentid = parent->pdg_id();}
0512         //cout << "      parent id = " << parentid << endl;
0513         HepMC::GenVertex* elevtx = (*p)->production_vertex();
0514         if(elevtx == jpsidecay_vtx) {
0515           cout << "electron from J/psi: " << pid << " " << status << " " << pt << " " << mass << " " << eta << endl;
0516           if(fabs(eta)<1.0) einacc = true;
0517           if(fabs(eta)<1.0 && pt>2.0) einacc2 = true;
0518           pt1 = pt;
0519           eta1 = eta;
0520         }
0521       }
0522       if(pid==11  && status==1 && pt>0.0) {
0523         HepMC::GenVertex* posvtx = (*p)->production_vertex();
0524         if(posvtx == jpsidecay_vtx) {
0525           cout << "positroni form J/psi: " << pid << " " << status << " " << pt << " " << mass << " " << eta << endl;
0526           if(fabs(eta)<1.0) pinacc = true;
0527           if(fabs(eta)<1.0 && pt>2.0) pinacc2 = true;
0528           pt2 = pt;
0529           eta2 = eta;
0530         }
0531       }
0532     }
0533 
0534     float type = 0.;
0535     if(einacc && pinacc) { type = 1.; }
0536     if(einacc2 && pinacc2) { type = 2.; }
0537 
0538 //  ntpmc2 = new TNtuple("ntpmc2","","type:mass:pt:eta:rap:pt1:pt2:eta1:eta2");
0539     float tmp[99];
0540     tmp[0] = type;
0541     tmp[1] = 0.;
0542     tmp[2] = jpsipt;
0543     tmp[3] = jpsieta;
0544     tmp[4] = jpsirap;
0545     tmp[5] = pt1;
0546     tmp[6] = pt2;
0547     tmp[7] = eta1;
0548     tmp[8] = eta2;
0549        ntpmc2->Fill(tmp);
0550 
0551 
0552   return Fun4AllReturnCodes::EVENT_OK;
0553 }
0554 
0555 //=====================================================================
0556 
0557 int sPHAnalysis::process_event_notracking(PHCompositeNode *topNode) {
0558 
0559   EventNumber++;
0560   int howoften = 1;
0561   if((EventNumber-1)%howoften==0) {
0562     cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
0563   }
0564   float tmp1[99];
0565 
0566   GlobalVertexMap *global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0567   //cout << "Number of GlobalVertexMap entries = " << global_vtxmap->size() << endl;
0568   double Zvtx = 0.;
0569   for (GlobalVertexMap::Iter iter = global_vtxmap->begin(); iter != global_vtxmap->end(); ++iter)
0570   {
0571     GlobalVertex *vtx = iter->second;
0572     if(vtx->get_id()==1) Zvtx = vtx->get_z(); // BBC vertex
0573     //cout << "Global vertex: " << vtx->get_id() << " " << vtx->get_z() << " " << vtx->get_t() << endl;
0574   }
0575   cout << "Global vertex Z = " << Zvtx << endl;
0576 
0577   RawClusterContainer* cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0578   if(!cemc_clusters) {
0579     cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
0580     return Fun4AllReturnCodes::ABORTEVENT;
0581   }
0582   cout << "   Number of CEMC clusters = " << cemc_clusters->size() << endl;
0583 
0584   RawTowerGeomContainer* _geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0585   if(!_geomEM) std::cerr<<"No TOWERGEOM_CEMC"<<std::endl;
0586   RawTowerGeomContainer* _geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0587   if(!_geomIH) std::cerr<<"No TOWERGEOM_HCALIN"<<std::endl;
0588   RawTowerGeomContainer* _geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0589   if(!_geomIH) std::cerr<<"No TOWERGEOM_HCALIN"<<std::endl;
0590 
0591   RawTowerContainer* _towersRawEM = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0592   if (!_towersRawEM) std::cerr<<"No TOWER_CALIB_CEMC Node"<<std::endl;
0593   RawTowerContainer* _towersRawIH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0594   if (!_towersRawIH) std::cerr<<"No TOWER_CALIB_HCALIN Node"<<std::endl;
0595   RawTowerContainer* _towersRawOH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0596   if (!_towersRawOH) std::cerr<<"No TOWER_CALIB_HCALOUT Node"<<std::endl;
0597 
0598   vector<TLorentzVector> electrons;
0599   vector<double> vhoe;
0600 
0601     // loop over CEMC clusters with E > 2 GeV 
0602     RawClusterContainer::Range begin_end = cemc_clusters->getClusters();
0603     RawClusterContainer::Iterator iter;
0604     for (iter = begin_end.first; iter != begin_end.second; ++iter)
0605     {
0606       RawCluster* cluster = iter->second;
0607 //      if(!cluster) { cout << "ERROR: bad cluster pointer = " << cluster << endl; continue; }
0608 //      else {
0609         double cemc_ecore = cluster->get_ecore();
0610         if(cemc_ecore<2.0) continue;
0611         double cemc_x = cluster->get_x();
0612         double cemc_y = cluster->get_y();
0613         double cemc_z = cluster->get_z() - Zvtx; // correct for event vertex position
0614         //double cemc_r = cluster->get_r();
0615         double cemc_r = sqrt(pow(cemc_x,2)+pow(cemc_y,2));
0616         double cemc_rr = sqrt(pow(cemc_r,2)+pow(cemc_z,2));
0617         double cemc_eta = asinh( cemc_z / cemc_r );
0618         double cemc_phi = atan2( cemc_y, cemc_x );
0619         double cemc_px = cemc_ecore * (cemc_x/cemc_rr);
0620         double cemc_py = cemc_ecore * (cemc_y/cemc_rr);
0621         double cemc_pz = cemc_ecore * (cemc_z/cemc_rr);
0622         //cout << "CEMC cluster: " << cemc_ecore << " " << cemc_eta << " " << cemc_phi << endl;
0623 
0624         // find closest CEMC tower in eta-phi space
0625           double distem = 99999.;
0626           RawTower* thetowerem = nullptr;
0627           RawTowerContainer::ConstRange begin_end_rawEM = _towersRawEM->getTowers();
0628           for (RawTowerContainer::ConstIterator rtiter = begin_end_rawEM.first; rtiter != begin_end_rawEM.second; ++rtiter) {
0629             RawTower *tower = rtiter->second;
0630             RawTowerGeom *tower_geom = _geomEM->get_tower_geometry(tower->get_key());
0631             double cemc_tower_phi  = tower_geom->get_phi();
0632             double cemc_tower_x = tower_geom->get_center_x();
0633             double cemc_tower_y = tower_geom->get_center_y();
0634             double cemc_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
0635             double cemc_tower_r = sqrt(pow(cemc_tower_x,2)+pow(cemc_tower_y,2));
0636             double cemc_tower_eta = asinh( cemc_tower_z / cemc_tower_r );
0637             double tmpdist = sqrt(pow(cemc_eta-cemc_tower_eta,2)+pow(cemc_phi-cemc_tower_phi,2));
0638             if(tmpdist<distem) { distem = tmpdist; thetowerem = tower; }
0639           }
0640           RawTowerGeom *thetower_geom_em = _geomEM->get_tower_geometry(thetowerem->get_key());
0641           unsigned int ietaem = thetower_geom_em->get_bineta();
0642           unsigned int jphiem = thetower_geom_em->get_binphi();
0643           h_notracking_etabins_em->Fill(double(ietaem));
0644           h_notracking_phibins_em->Fill(double(jphiem));
0645 
0646 
0647 
0648           // find closest HCALIN tower in eta-phi space
0649           double distin = 99999.;
0650           RawTower* thetowerin = nullptr;
0651           RawTowerContainer::ConstRange begin_end_rawIN = _towersRawIH->getTowers();
0652           for (RawTowerContainer::ConstIterator rtiter = begin_end_rawIN.first; rtiter != begin_end_rawIN.second; ++rtiter) {
0653             RawTower *tower = rtiter->second;
0654             RawTowerGeom *tower_geom = _geomIH->get_tower_geometry(tower->get_key());
0655             double hcalin_tower_phi  = tower_geom->get_phi();
0656             double hcalin_tower_x = tower_geom->get_center_x();
0657             double hcalin_tower_y = tower_geom->get_center_y();
0658             double hcalin_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
0659             double hcalin_tower_r = sqrt(pow(hcalin_tower_x,2)+pow(hcalin_tower_y,2));
0660             double hcalin_tower_eta = asinh( hcalin_tower_z / hcalin_tower_r );
0661             double tmpdist = sqrt(pow(cemc_eta-hcalin_tower_eta,2)+pow(cemc_phi-hcalin_tower_phi,2));
0662             if(tmpdist<distin) { distin = tmpdist; thetowerin = tower; }
0663           }
0664           RawTowerGeom *thetower_geom = _geomIH->get_tower_geometry(thetowerin->get_key());
0665           unsigned int ieta = thetower_geom->get_bineta();
0666           unsigned int jphi = thetower_geom->get_binphi();
0667           h_notracking_etabins->Fill(double(ieta));
0668           h_notracking_phibins->Fill(double(jphi));
0669 
0670           // Calcuate 3x3 energy deposit in HCALIN. Inner HCAL has 24 bins in eta and 64 bins in phi
0671           double e3x3in = 0.;
0672           if(ieta<1 || ieta>22) continue; // ignore clusters in edge towers
0673           for(unsigned int i=0; i<=2; i++) {
0674             for(unsigned int j=0; j<=2; j++) {
0675               unsigned int itmp = ieta-1+i;
0676               unsigned int jtmp = 0;  
0677               if(jphi==0 && j==0) { jtmp = 63; } // wrap around
0678               else if(jphi==63 && j==2) { jtmp = 0; } // wrap around
0679               else { jtmp = jphi-1+j; } 
0680               RawTower* tmptower = _towersRawIH->getTower(itmp,jtmp);
0681               if(tmptower) { e3x3in += tmptower->get_energy(); }
0682             }
0683           }
0684           h_notracking_eoh->Fill(e3x3in/cemc_ecore);
0685 
0686           // find closest HCALOUT tower in eta-phi space
0687           double distout = 99999.;
0688           RawTower* thetowerout = nullptr;
0689           RawTowerContainer::ConstRange begin_end_rawON = _towersRawOH->getTowers();
0690           for (RawTowerContainer::ConstIterator rtiter = begin_end_rawON.first; rtiter != begin_end_rawON.second; ++rtiter) {
0691             RawTower *tower = rtiter->second;
0692             RawTowerGeom *tower_geom = _geomOH->get_tower_geometry(tower->get_key());
0693             double hcalout_tower_phi  = tower_geom->get_phi();
0694             double hcalout_tower_x = tower_geom->get_center_x();
0695             double hcalout_tower_y = tower_geom->get_center_y();
0696             double hcalout_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
0697             double hcalout_tower_r = sqrt(pow(hcalout_tower_x,2)+pow(hcalout_tower_y,2));
0698             double hcalout_tower_eta = asinh( hcalout_tower_z / hcalout_tower_r );
0699             double tmpdist = sqrt(pow(cemc_eta-hcalout_tower_eta,2)+pow(cemc_phi-hcalout_tower_phi,2));
0700             if(tmpdist<distout) { distout = tmpdist; thetowerout = tower; }
0701           }
0702           RawTowerGeom *thetower_geomout = _geomOH->get_tower_geometry(thetowerout->get_key());
0703           unsigned int ietaout = thetower_geomout->get_bineta();
0704           unsigned int jphiout = thetower_geomout->get_binphi();
0705           h_notracking_etabinsout->Fill(double(ietaout));
0706           h_notracking_phibinsout->Fill(double(jphiout));
0707 
0708           double e3x3out = 0.;
0709           if(ietaout<1 || ietaout>22) continue; // ignore clusters in edge towers
0710           for(unsigned int i=0; i<=2; i++) {
0711             for(unsigned int j=0; j<=2; j++) {
0712               unsigned int itmp = ietaout-1+i;
0713               unsigned int jtmp = 0;
0714               if(jphiout==0 && j==0) { jtmp = 63; } // wrap around
0715               else if(jphiout==63 && j==2) { jtmp = 0; } // wrap around
0716               else { jtmp = jphiout-1+j; }
0717               RawTower* tmptower = _towersRawOH->getTower(itmp,jtmp);
0718               if(tmptower) { e3x3out += tmptower->get_energy(); }
0719             }
0720           }
0721 
0722 
0723           if(e3x3in/cemc_ecore>0.06) continue; // reject hadrons, 90% eID efficiency is with 0.058 cut
0724 
0725           TLorentzVector tmp = TLorentzVector(cemc_px,cemc_py,cemc_pz,cemc_ecore);
0726           electrons.push_back(tmp);
0727           vhoe.push_back(e3x3in/cemc_ecore);
0728 
0729 //      } // valid CEMC cluster
0730     } // end loop over CEMC clusters
0731     cout << "number of selected electrons = " << electrons.size() << " " << vhoe.size() << endl;
0732 
0733 // Make Upsilons
0734 
0735   if(electrons.size()>=2) {
0736     for(long unsigned int i=0; i<electrons.size()-1; i++) {
0737       for(long unsigned int j=i+1; j<electrons.size(); j++) {
0738         TLorentzVector pair = electrons[i]+electrons[j];
0739         cout << i << " " << j << endl;
0740         tmp1[0] = pair.M();
0741         tmp1[1] = pair.Pt();
0742         cout << pair.M() << " " << pair.Pt() << endl;
0743         tmp1[2] = (electrons[i]).Pt();
0744         tmp1[3] = (electrons[j]).Pt();
0745         cout << vhoe[i] << " " << vhoe[j] << endl;
0746         tmp1[4] = vhoe[i];
0747         tmp1[5] = vhoe[j];
0748         cout << "filling..." << endl;
0749           ntp_notracking->Fill(tmp1);
0750         cout << "done." << endl;
0751       }
0752     }
0753   }
0754 
0755   cout << "returning..." << endl;
0756   return Fun4AllReturnCodes::EVENT_OK;
0757 
0758 }
0759 
0760 //======================================================================
0761 
0762 double sPHAnalysis::Get_CAL_e3x3(SvtxTrack* track, RawTowerContainer* _towersRawOH, RawTowerGeomContainer* _geomOH, int what, double Zvtx, double &dphi, double &deta) {
0763 
0764 double e3x3 = 0.;
0765 double pathlength = 999.;
0766 vector<double> proj;
0767 for (SvtxTrack::StateIter stateiter = track->begin_states(); stateiter != track->end_states(); ++stateiter)
0768 {
0769   SvtxTrackState *trackstate = stateiter->second;
0770   if(trackstate) { proj.push_back(trackstate->get_pathlength()); }
0771 }
0772 if(what==0)      { pathlength = proj[proj.size()-3]; } // CEMC
0773 else if(what==1) { pathlength = proj[proj.size()-2]; } // HCALIN
0774 else if(what==2) { pathlength = proj[proj.size()-1]; } // HCALOUT
0775 else { dphi = 9999.; deta = 9999.; return e3x3;}
0776 
0777   double track_eta = 999.;
0778   double track_phi = 999.;
0779   SvtxTrackState* trackstate = track->get_state(pathlength);
0780   if(trackstate) {
0781     double track_x = trackstate->get_x();
0782     double track_y = trackstate->get_y();
0783     double track_z = trackstate->get_z() - Zvtx;
0784     double track_r = sqrt(track_x*track_x+track_y*track_y);
0785     track_eta = asinh( track_z / track_r );
0786     track_phi = atan2( track_y, track_x );
0787   } else { cout << "track state not found!" << endl; dphi = 9999.; deta = 9999.; return e3x3; }
0788 
0789   double dist = 9999.;
0790   RawTower* thetower = nullptr;
0791   RawTowerContainer::ConstRange begin_end_rawOH = _towersRawOH->getTowers();
0792   for (RawTowerContainer::ConstIterator rtiter = begin_end_rawOH.first; rtiter != begin_end_rawOH.second; ++rtiter) {
0793     RawTower *tower = rtiter->second;
0794     RawTowerGeom *tower_geom = _geomOH->get_tower_geometry(tower->get_key());
0795     //double tower_phi  = tower_geom->get_phi();
0796     double tower_x = tower_geom->get_center_x();
0797     double tower_y = tower_geom->get_center_y();
0798     double tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
0799     double tower_r = sqrt(pow(tower_x,2)+pow(tower_y,2));
0800     double tower_eta = asinh( tower_z / tower_r );
0801     double tower_phi = atan2( tower_y , tower_x );
0802     double tmpdist = sqrt(pow(track_eta-tower_eta,2)+pow(track_phi-tower_phi,2));
0803     if(tmpdist<dist) { dist = tmpdist; thetower = tower; deta = fabs(track_eta-tower_eta); dphi = fabs(track_phi-tower_phi); }
0804   }
0805   cout << "dist: " << dist << " " << deta << " " << dphi << endl;
0806 
0807   if(!thetower) { dphi = 9999.; deta = 9999.; return e3x3; }
0808   RawTowerGeom *thetower_geom = _geomOH->get_tower_geometry(thetower->get_key());
0809   unsigned int ieta = thetower_geom->get_bineta();
0810   unsigned int jphi = thetower_geom->get_binphi();
0811 
0812   unsigned int maxbinphi = 63; if(what==0) maxbinphi = 255;
0813   unsigned int maxbineta = 23; if(what==0) maxbineta = 93;
0814     // calculate e3x3
0815     for(unsigned int i=0; i<=2; i++) {
0816       for(unsigned int j=0; j<=2; j++) {
0817         unsigned int itmp = ieta-1+i;
0818         unsigned int jtmp = 0;
0819         if(jphi==0 && j==0) { jtmp = maxbinphi; }      // wrap around
0820         else if(jphi==maxbinphi && j==2) { jtmp = 0; } // wrap around
0821         else { jtmp = jphi-1+j; }
0822         if(itmp>=0 && itmp<=maxbineta) { 
0823           RawTower* tmptower = _towersRawOH->getTower(itmp,jtmp);
0824           if(tmptower) { e3x3 += tmptower->get_energy(); }
0825         }
0826       }
0827     }
0828 
0829   return e3x3;
0830 }
0831 
0832 //======================================================================
0833 
0834 int sPHAnalysis::process_event_filtered(PHCompositeNode *topNode) {
0835 EventNumber++;
0836 
0837   cout << "Event # " << EventNumber << endl;
0838 
0839   SvtxTrackMap* _trackmap = nullptr;
0840   TrackSeedContainer* _trackseedcontainer_svtx = nullptr;
0841   TrackSeedContainer* _trackseedcontainer_silicon = nullptr;
0842   TrackSeedContainer* _trackseedcontainer_tpc = nullptr;
0843   TrkrClusterContainer* _trkrclusters = nullptr;
0844   RawClusterContainer* _cemc_clusters = nullptr;
0845 
0846   _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0847   if(!_trackmap) { cerr << PHWHERE << "ERROR: SvtxTrackMap node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0848   cout << "   Number of tracks = " << _trackmap->size() << endl;
0849 
0850   _trkrclusters  = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0851   if(!_trkrclusters) { cerr << PHWHERE << "ERROR: TRKR_CLUSTER node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0852   cout << "   Number of TRKR clusters = " << _trkrclusters->size() << endl;
0853 
0854   _trackseedcontainer_svtx = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
0855   if(!_trackseedcontainer_svtx) { cerr << PHWHERE << "ERROR: SvtxTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0856 
0857   _trackseedcontainer_silicon = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0858   if(!_trackseedcontainer_silicon) { cerr << PHWHERE << "ERROR: SiliconTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0859 
0860   _trackseedcontainer_tpc = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0861   if(!_trackseedcontainer_tpc) { cerr << PHWHERE << "ERROR: TpcTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0862 
0863   _cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0864   if(!_cemc_clusters) { cerr << PHWHERE << "ERROR: CLUSTER_CEMC node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0865   cout << "   Number of CEMC clusters = " << _cemc_clusters->size() << endl;
0866 
0867   return 0;
0868 }
0869 
0870 
0871 
0872 //======================================================================
0873 
0874 int sPHAnalysis::process_event_test(PHCompositeNode *topNode) {
0875 
0876   int evtno = EventNumber;
0877   int evtno5 = evtno%5;
0878 
0879   EventNumber++;
0880   int howoften = 1; 
0881   if((EventNumber-1)%howoften==0) { 
0882     cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
0883   }
0884   //if(EventNumber==1) topNode->print();
0885   float tmp1[99];
0886   float tmpb[99];
0887 
0888   GlobalVertexMap *global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0889   //cout << "Number of GlobalVertexMap entries = " << global_vtxmap->size() << endl;
0890   double Zvtx = 0.;
0891   for (GlobalVertexMap::Iter iter = global_vtxmap->begin(); iter != global_vtxmap->end(); ++iter)
0892   {
0893     GlobalVertex *vtx = iter->second;
0894     if(vtx->get_id()==1) Zvtx = vtx->get_z(); // BBC vertex
0895     //cout << "Global vertex: " << vtx->get_id() << " " << vtx->get_z() << " " << vtx->get_t() << endl;
0896   }
0897   cout << "Global BBC vertex Z = " << Zvtx << endl;
0898 
0899 // TRUTH ------------------------------------------------------------------------
0900   float impact_parameter = 999.;
0901   int npart = 0;
0902   int ginacc = 0;
0903 /*
0904 
0905   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0906   PHHepMCGenEvent *genevt = nullptr;
0907   if(!genevtmap) {cout << "NO PHHepMCGenEventMap node found!" << endl;
0908 
0909   for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
0910   {
0911     genevt = iter->second;
0912     if(genevt->get_embedding_id()==0) break;
0913     if(!genevt) {cout<<"ERROR: no PHHepMCGenEvent!" << endl; return Fun4AllReturnCodes::ABORTEVENT;}
0914   }
0915 
0916   HepMC::GenEvent *event = genevt->getEvent();
0917   if (!event) { cout << PHWHERE << "ERROR: no HepMC::GenEvent!" << endl; return Fun4AllReturnCodes::ABORTEVENT;}
0918   npart = event->particles_size();
0919   int nvert = event->vertices_size();
0920   cout << "HepMC::GenEvent: Number of particles, vertuces = " << npart << " " << nvert << endl;
0921 
0922   HepMC::HeavyIon* hi = event->heavy_ion();
0923   if(hi) { 
0924     impact_parameter = hi->impact_parameter(); 
0925     cout << "HepMC::GenEvent: impact parameter = " << impact_parameter << endl;
0926   }
0927 
0928   PHG4TruthInfoContainer* truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0929   if(!truth_container) {
0930     cerr << PHWHERE << " ERROR: Can not find G4TruthInfo node." << endl;
0931     return Fun4AllReturnCodes::ABORTEVENT;
0932   }
0933   PHG4TruthInfoContainer::ConstRange range = truth_container->GetPrimaryParticleRange();
0934 
0935   vector<TLorentzVector> gparticles;
0936 
0937 
0938     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0939     {
0940       PHG4Particle* g4particle = iter->second;
0941       int gflavor  = g4particle->get_pid();
0942       double gmass = 0.;
0943       if(fabs(gflavor)==11)       { gmass = 0.000511; } 
0944       else if(fabs(gflavor)==211) { gmass = 0.13957; } 
0945       else if(fabs(gflavor)==321) { gmass = 0.49368; } 
0946       else if(fabs(gflavor)==2212) { gmass = 0.93827; } 
0947       else { continue; }
0948         double gpx = g4particle->get_px();
0949         double gpy = g4particle->get_py();
0950         double gpz= g4particle->get_pz();
0951         double gpt = sqrt(gpx*gpx+gpy*gpy);
0952         double ge = sqrt(gpt*gpt + gpz*gpz + gmass * gmass);
0953         TLorentzVector tmp = TLorentzVector(gpx,gpy,gpz,ge);
0954         if(gpt>0.5 && fabs(tmp.Eta())<1.0) { ginacc++; }
0955 
0956       int trackid = g4particle->get_track_id();
0957       if(trackid>truth_container->GetNumPrimaryVertexParticles()-50) { //embedded particles are the last ones
0958       //if(trackid>truth_container->GetNumPrimaryVertexParticles()-2) { // Quarkonia
0959         double gpx = g4particle->get_px();
0960         double gpy = g4particle->get_py();
0961         double gpz= g4particle->get_pz();
0962         double gpt = sqrt(gpx*gpx+gpy*gpy);
0963         //double phi = atan2(gpy,gpx);
0964         //double eta = asinh(gpz/gpt);
0965         double ge = sqrt(gpt*gpt + gpz*gpz + gmass * gmass);
0966         //int primid =  g4particle->get_primary_id();
0967         //int parentid = g4particle->get_parent_id();
0968         TLorentzVector tmp = TLorentzVector(gpx,gpy,gpz,ge);
0969         gparticles.push_back(tmp);
0970         //cout << "embedded: " << gflavor << " " << gpt << " " << gmass << endl;
0971       }
0972     }
0973     cout << "number of embedded particles = " << gparticles.size() << endl;
0974 
0975 //    TLorentzVector ele = gparticles[0];
0976 //    TLorentzVector pos = gparticles[1];
0977 //    TLorentzVector jpsi = ele + pos;;
0978 //    cout << "Embedded J/psi: " << jpsi.Pt() << " " << jpsi.Eta() << endl; 
0979 */
0980 
0981 // RECO -------------------------------------------------------------------
0982 
0983   SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0984   if(!trackmap) {
0985     cerr << PHWHERE << " ERROR: Can not find SvtxTrackMap node." << endl;
0986     return Fun4AllReturnCodes::ABORTEVENT;
0987   }
0988   cout << "   Number of tracks = " << trackmap->size() << endl;
0989   double mult = trackmap->size();
0990   hmult->Fill(mult);
0991 
0992   RawClusterContainer* cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0993   if(!cemc_clusters) {
0994     cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
0995     return Fun4AllReturnCodes::ABORTEVENT;
0996   }
0997   cout << "   Number of CEMC clusters = " << cemc_clusters->size() << endl;
0998   int cemcmult = cemc_clusters->size();
0999 
1000   RawClusterContainer* hcalin_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALIN");
1001   if(!hcalin_clusters) {
1002     cerr << PHWHERE << " ERROR: Can not find CLUSTER_HCALIN node." << endl;
1003     return Fun4AllReturnCodes::ABORTEVENT;
1004   }
1005   RawClusterContainer* hcalout_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_HCALOUT");
1006   if(!hcalout_clusters) {
1007     cerr << PHWHERE << " ERROR: Can not find CLUSTER_HCALOUT node." << endl;
1008     return Fun4AllReturnCodes::ABORTEVENT;
1009   }
1010 
1011   RawTowerGeomContainer* _geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
1012   if(!_geomEM) std::cerr<<"No TOWERGEOM_CEMC"<<std::endl;
1013   RawTowerGeomContainer* _geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
1014   if(!_geomIH) std::cerr<<"No TOWERGEOM_HCALIN"<<std::endl;
1015   RawTowerGeomContainer* _geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
1016   if(!_geomOH) std::cerr<<"No TOWERGEOM_HCALIN"<<std::endl;
1017 
1018   RawTowerContainer* _towersEM = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
1019   if (!_towersEM) std::cerr<<"No TOWER_CALIB_CEMC Node"<<std::endl;
1020   RawTowerContainer* _towersIH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
1021   if (!_towersIH) std::cerr<<"No TOWER_CALIB_HCALIN Node"<<std::endl;
1022   RawTowerContainer* _towersOH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
1023   if (!_towersOH) std::cerr<<"No TOWER_CALIB_HCALOUT Node"<<std::endl;
1024   
1025 //------------------------------------------------------------------------------------------
1026 /*
1027     RawClusterContainer::Range begin_end = cemc_clusters->getClusters();
1028     RawClusterContainer::Iterator iter;
1029     for (iter = begin_end.first; iter != begin_end.second; ++iter)
1030     {
1031       RawCluster* cluster = iter->second;
1032       if(!cluster) { cout << "ERROR: bad cluster pointer = " << cluster << endl; continue; }
1033       else {
1034         double cemc_ecore = cluster->get_ecore();
1035         //double cemc_e = cluster->get_energy();
1036         //double cemc_phi = cluster->get_phi();
1037         double cemc_x = cluster->get_x();
1038         double cemc_y = cluster->get_y();
1039         double cemc_z = cluster->get_z();
1040         double cemc_r = cluster->get_r();
1041         double cemc_eta = asinh( cemc_z / cemc_r );
1042         double cemc_phi = atan2( cemc_y, cemc_x );
1043       }
1044     }
1045 */
1046 
1047 //------------------------------------------------------------------------------------------
1048 
1049 int ngood = 0;
1050 cout << "starting loop over tracks..." << endl;
1051   for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
1052   {
1053     SvtxTrack *track = iter->second;
1054     if(!track) { cout << "ERROR: bad track pointer = " << track << endl; continue; }
1055 
1056     double charge = track->get_charge();
1057     double px = track->get_px();
1058     double py = track->get_py();
1059     double pz = track->get_pz();
1060     double pt = sqrt(px * px + py * py);
1061     double mom = sqrt(px * px + py * py + pz * pz);
1062     double eta = track->get_eta();
1063     //double phi = track->get_phi();
1064     double chisq = track->get_chisq();
1065     double ndf = track->get_ndf();
1066     double chi2 = 999.;
1067       if(ndf!=0.) chi2 = chisq/ndf;
1068     double dca3d_xy = track->get_dca3d_xy();
1069     if(pt>0.5 && chi2<5.) { ngood++; }
1070       if(pt<2.0) continue;
1071 
1072 // Find matching CEMC cluster
1073     double cemc_dphi = 99999.;
1074     double cemc_deta = 99999.;
1075     RawCluster* clus = MatchClusterCEMC(track,cemc_clusters, cemc_dphi, cemc_deta, Zvtx, 1);
1076     double cemc_ecore = 0.;
1077     double cemc_prob = 0.;
1078     double cemc_chi2 = 0.;
1079     if(clus) {
1080       cemc_ecore = clus->get_ecore();
1081       cemc_prob = clus->get_prob();
1082     }
1083 
1084 // Calculate e3x3 for calorimeters
1085     double cemc_deta2 = 9999.;
1086     double cemc_dphi2 = 9999.;
1087     double cemc_e3x3 = Get_CAL_e3x3(track, _towersEM, _geomEM, 0, Zvtx, cemc_dphi2, cemc_deta2);
1088     double hcalin_deta = 9999.;
1089     double hcalin_dphi = 9999.;
1090     double hcalin_e3x3 = Get_CAL_e3x3(track, _towersIH, _geomIH, 1, Zvtx, hcalin_dphi, hcalin_deta);
1091     double hcalout_deta = 9999.;
1092     double hcalout_dphi = 9999.;
1093     double hcalout_e3x3 = Get_CAL_e3x3(track, _towersOH, _geomOH, 2, Zvtx, hcalout_dphi, hcalout_deta);
1094 
1095     double hcalin_clusdphi = 9999.;
1096     double hcalin_clusdeta = 9999.;
1097     double hcalout_clusdphi = 9999.;
1098     double hcalout_clusdeta = 9999.;
1099     RawCluster* clus_hcalin = MatchClusterCEMC(track, hcalin_clusters, hcalin_clusdphi, hcalin_clusdeta, Zvtx, 2);
1100     RawCluster* clus_hcalout = MatchClusterCEMC(track, hcalout_clusters, hcalout_clusdphi, hcalout_clusdeta, Zvtx, 3);
1101     double hcalin_cluse = clus_hcalin->get_energy();
1102     double hcalout_cluse = clus_hcalout->get_energy();
1103 
1104     cout << "track: " << charge << " " << pt << "    " << cemc_ecore << " " << cemc_e3x3 << " " << hcalin_e3x3 << " " << hcalout_e3x3 << endl;
1105 
1106 // count hits
1107     int nmvtx = 0; int nintt = 0; int ntpc = 0;
1108     auto siseed = track->get_silicon_seed();
1109     if(siseed) {
1110       for (auto iter = siseed->begin_cluster_keys(); iter != siseed->end_cluster_keys(); ++iter) {
1111         TrkrDefs::cluskey cluster_key = *iter;
1112         auto trkrid = TrkrDefs::getTrkrId(cluster_key);
1113         if(trkrid==TrkrDefs::mvtxId) nmvtx++;
1114         if(trkrid==TrkrDefs::inttId) nintt++;
1115         //int layer = TrkrDefs::getLayer(cluster_key);
1116         //if(0<=layer && layer<=2) nmvtx++;
1117         //if(3<=layer && layer<=6) nintt++;
1118       }
1119     }
1120     auto tpcseed = track->get_tpc_seed();
1121     if(tpcseed) {
1122       for (auto iter = tpcseed->begin_cluster_keys(); iter != tpcseed->end_cluster_keys(); ++iter) {
1123         TrkrDefs::cluskey cluster_key = *iter;
1124         int layer = TrkrDefs::getLayer(cluster_key);
1125         if(layer>6) ntpc++;
1126       }
1127     }
1128     cout << "    ntpc, nmvtx, nitt= " << ntpc << " " << nmvtx << " " << nintt << endl;
1129       
1130 // print all track states
1131 /*
1132   vector<double> proj;
1133   for (SvtxTrack::StateIter stateiter = track->begin_states(); stateiter != track->end_states(); ++stateiter)
1134   {
1135     SvtxTrackState *trackstate = stateiter->second;
1136     if(trackstate) { 
1137       cout << "state: " << trackstate->get_pathlength() << endl; 
1138       proj.push_back(trackstate->get_pathlength());
1139     } 
1140   }
1141   //cout << "projection radii = " << proj[proj.size()-3] << " " << proj[proj.size()-2] << " " << proj[proj.size()-1] << endl;
1142 */
1143     double gdist = 999.;
1144 /*    
1145     for(unsigned int i=0; i<gparticles.size(); i++) {
1146       double gphi = gparticles[i].Phi();
1147       double geta = gparticles[i].Eta();
1148       double tmpdist = sqrt(pow(phi-gphi,2)+pow(eta-geta,2));
1149       if(tmpdist<gdist) { gdist = tmpdist; }
1150     }
1151     cout << "gdist: " << gdist << endl;
1152 //    if(gdist>0.001) continue;
1153 */
1154 
1155     tmp1[0] = charge;
1156     tmp1[1] = pt;
1157     tmp1[2] = eta;
1158     tmp1[3] = mom;
1159     tmp1[4] = nmvtx;
1160     tmp1[5] = ntpc;
1161     tmp1[6] = chisq;
1162     tmp1[7] = ndf;
1163     tmp1[8] = cemc_ecore;
1164     tmp1[9] = cemc_e3x3;
1165     tmp1[10] = cemc_prob;
1166     tmp1[11] = cemc_chi2;
1167     tmp1[12] = cemc_dphi;
1168     tmp1[13] = cemc_deta;
1169     tmp1[14] = chi2;
1170     tmp1[15] = dca3d_xy;
1171     tmp1[16] = hcalin_e3x3;
1172     tmp1[17] = hcalout_e3x3;
1173     tmp1[18] = gdist;
1174     tmp1[19] = cemc_dphi2;
1175     tmp1[20] = cemc_dphi2;
1176     tmp1[21] = hcalin_dphi;
1177     tmp1[22] = hcalin_deta;
1178     tmp1[23] = hcalout_deta;
1179     tmp1[24] = hcalout_deta;
1180     tmp1[25] = hcalin_clusdphi;
1181     tmp1[26] = hcalin_clusdeta;
1182     tmp1[27] = hcalin_cluse;
1183     tmp1[28] = hcalout_clusdphi;
1184     tmp1[29] = hcalout_clusdeta;
1185     tmp1[30] = hcalout_cluse;
1186       ntppid->Fill(tmp1);
1187 
1188   } // end loop over tracks
1189 
1190 //  ntpb = new TNtuple("ntpb","","bimp:mult:truemult:cemcmult");
1191   tmpb[0] = impact_parameter;
1192   tmpb[1] = mult;
1193   tmpb[2] = npart;
1194   tmpb[3] = cemcmult;
1195   tmpb[4] = float(evtno);
1196   tmpb[5] = float(evtno5);
1197   tmpb[6] = float(ginacc);
1198   tmpb[7] = float(ngood);
1199     ntpb->Fill(tmpb);
1200 
1201   return Fun4AllReturnCodes::EVENT_OK;
1202 }
1203 //================================================================================
1204 
1205 TVector3 sPHAnalysis::GetProjectionToCalorimeter(SvtxTrack* track, int what) {
1206 // what = 1 CEMC
1207 // what = 2 HCALIN
1208 // what = 3 HCALOUT
1209 
1210   TVector3 projection; // 0,0,0
1211 
1212   vector<double> proj;
1213   for (SvtxTrack::StateIter stateiter = track->begin_states(); stateiter != track->end_states(); ++stateiter)
1214   {
1215     SvtxTrackState *trackstate = stateiter->second;
1216     if(trackstate) { proj.push_back(trackstate->get_pathlength()); }
1217   }
1218   double pathlength = proj[proj.size()+4-what]; // CEMC is next next to last, usually 93.5
1219 
1220   SvtxTrackState* trackstate = track->get_state(pathlength); // at CEMC inner face
1221   if(trackstate) {
1222     double track_x = trackstate->get_x();
1223     double track_y = trackstate->get_y();
1224     double track_z = trackstate->get_z();
1225     projection.SetX(track_x);
1226     projection.SetY(track_y);
1227     projection.SetZ(track_z);
1228   }
1229 
1230   return projection;
1231 }
1232 
1233 //=================================================================================
1234 
1235 RawCluster* sPHAnalysis::MatchClusterCEMC(SvtxTrack* track, RawClusterContainer* cemc_clusters, double &dphi, double &deta, double Zvtx, int what) {
1236 
1237   RawCluster* returnCluster = NULL;
1238   double track_eta = 99999.;
1239   double track_phi = 99999.;
1240   dphi = 99999.;
1241   deta = 99999.;
1242 
1243   vector<double> proj;
1244   for (SvtxTrack::StateIter stateiter = track->begin_states(); stateiter != track->end_states(); ++stateiter)
1245   {
1246     SvtxTrackState *trackstate = stateiter->second;
1247     if(trackstate) { proj.push_back(trackstate->get_pathlength()); }
1248   }
1249   double pathlength = 0.;
1250   if(what==1) {
1251     pathlength = proj[proj.size()-3]; // CEMC is next next to last
1252   } else if(what==2) {
1253     pathlength = proj[proj.size()-2]; // HCALIN is next to last
1254   } else if(what==3) {
1255     pathlength = proj[proj.size()-1]; // HCALOUT is the last
1256   } else {return returnCluster; }
1257 
1258   SvtxTrackState* trackstate = track->get_state(pathlength); // at CEMC inner face
1259   if(trackstate) {
1260     double track_x = trackstate->get_x();
1261     double track_y = trackstate->get_y();
1262     double track_z = trackstate->get_z() - Zvtx;
1263     double track_r = sqrt(track_x*track_x+track_y*track_y);
1264     track_eta = asinh( track_z / track_r );
1265     track_phi = atan2( track_y, track_x );
1266   } else { return returnCluster; }
1267 
1268   if(track_eta == 99999. || track_phi == 99999.) { return returnCluster; }
1269   double dist = 99999.;
1270 
1271     RawClusterContainer::Range begin_end = cemc_clusters->getClusters();
1272     RawClusterContainer::Iterator iter;
1273     for (iter = begin_end.first; iter != begin_end.second; ++iter)
1274     {
1275       RawCluster* cluster = iter->second;
1276       if(!cluster) continue;
1277       else {
1278         double cemc_ecore = cluster->get_ecore();
1279         if(cemc_ecore<0.0) continue;
1280         double cemc_x = cluster->get_x();
1281         double cemc_y = cluster->get_y();
1282         double cemc_z = cluster->get_z() - Zvtx;
1283         double cemc_r = cluster->get_r();
1284         double cemc_eta = asinh( cemc_z / cemc_r );
1285         double cemc_phi = atan2( cemc_y, cemc_x );
1286         double tmpdist = sqrt(pow((cemc_eta-track_eta),2)+pow((cemc_phi-track_phi),2));
1287         if(tmpdist<dist) {
1288           dist = tmpdist; returnCluster = cluster; dphi = fabs(cemc_phi-track_phi); deta = fabs(cemc_eta-track_eta);
1289         }
1290       }
1291     }
1292 
1293   return returnCluster;
1294 }
1295 
1296 
1297 //=================================================================================
1298 
1299 int sPHAnalysis::process_event_upsilons(PHCompositeNode *topNode) {
1300   EventNumber++;
1301   int howoften = 1;
1302   if((EventNumber-1)%howoften==0) {
1303     cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
1304   }
1305   if(EventNumber==1) topNode->print();
1306 
1307 
1308   SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
1309   if(!trackmap) {
1310     cerr << PHWHERE << " ERROR: Can not find SvtxTrackMap node." << endl;
1311     return Fun4AllReturnCodes::ABORTEVENT;
1312   }
1313   cout << "   Number of tracks = " << trackmap->size() << endl;
1314 
1315   SvtxVertexMap *vtxmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
1316   if(!vtxmap) {
1317       cout << "SvtxVertexMap node not found!" << endl;
1318       return Fun4AllReturnCodes::ABORTEVENT;
1319   }
1320   cout << "Number of SvtxVertexMap entries = " << vtxmap->size() << endl;
1321 
1322   double Zvtx = 0.;
1323   GlobalVertexMap *global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
1324   if(!global_vtxmap) {
1325     cerr << PHWHERE << " ERROR: Can not find GlobalVertexMap node." << endl;
1326     return Fun4AllReturnCodes::ABORTEVENT;
1327   }
1328   for (GlobalVertexMap::Iter iter = global_vtxmap->begin(); iter != global_vtxmap->end(); ++iter)
1329   { GlobalVertex *vtx = iter->second; if(vtx->get_id()==1) { Zvtx = vtx->get_z(); } }
1330   cout << "Global BBC vertex Z = " << Zvtx << endl;
1331 
1332   RawClusterContainer* cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
1333   if(!cemc_clusters) {
1334     cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
1335     return Fun4AllReturnCodes::ABORTEVENT;
1336   }
1337   else { cout << "FOUND CLUSTER_CEMC node." << endl; }
1338 
1339   vector<TLorentzVector> electrons;
1340   vector<TLorentzVector> positrons;
1341   vector<double> vpchi2;
1342   vector<double> vmchi2;
1343   vector<double> vpdca;
1344   vector<double> vmdca;
1345   vector<double> vpeop3x3;
1346   vector<double> vmeop3x3;
1347   vector<int> vpnmvtx;
1348   vector<int> vpnintt;
1349   vector<int> vpntpc;
1350   vector<int> vmnmvtx;
1351   vector<int> vmnintt;
1352   vector<int> vmntpc;
1353 
1354   for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
1355   {
1356     SvtxTrack *track = iter->second;
1357     if(!track) {
1358       cout << "ERROR: bad track pointer = " << track << endl;
1359       continue;
1360     }
1361 
1362     double charge = track->get_charge();
1363     double px = track->get_px();
1364     double py = track->get_py();
1365     double pz = track->get_pz();
1366     double pt = sqrt(px * px + py * py);
1367     double mom = sqrt(pt*pt+pz*pz);
1368     double ee = sqrt(mom*mom+0.000511*0.000511);
1369       //if(pt<0.5) continue;
1370       if(pt<2.0) continue;
1371     //double phi = track->get_phi();
1372     double eta = track->get_eta();
1373     double chisq = track->get_chisq();
1374     double ndf = track->get_ndf();
1375     double chi2 = 999.;
1376     double dca3d_xy = track->get_dca3d_xy(); if(ndf!=0.) chi2 = chisq/ndf;
1377 
1378     double cemc_dphi = 99999.;
1379     double cemc_deta = 99999.;
1380     RawCluster* clus = MatchClusterCEMC(track, cemc_clusters, cemc_dphi, cemc_deta, Zvtx, 1);
1381     double cemc_ecore = 0.;
1382     if(clus) cemc_ecore = clus->get_ecore();
1383     if(cemc_ecore/mom<0.7) continue; // not an electron
1384     double trk_x = track->get_x();
1385     double trk_y = track->get_y();
1386     double trk_z = track->get_z();
1387     cout << "track: " << charge << " " << pt << " " << cemc_ecore/mom << " " << trk_x << " " << trk_y << " " << trk_z << endl;
1388 
1389     int nmvtx = 0; int nintt = 0; int ntpc = 0;
1390     auto siseed = track->get_silicon_seed();
1391     for (auto iter = siseed->begin_cluster_keys(); iter != siseed->end_cluster_keys(); ++iter) {
1392       TrkrDefs::cluskey cluster_key = *iter;
1393       //int layer = TrkrDefs::getLayer(cluster_key);
1394       auto trkrid = TrkrDefs::getTrkrId(cluster_key);
1395       if(trkrid==TrkrDefs::mvtxId) nmvtx++;
1396       if(trkrid==TrkrDefs::inttId) nintt++;
1397       //if(0<=layer && layer<=2) nmvtx++;
1398       //if(3<=layer && layer<=6) nintt++;
1399     }
1400     auto tpcseed = track->get_tpc_seed();
1401     for (auto iter = tpcseed->begin_cluster_keys(); iter != tpcseed->end_cluster_keys(); ++iter) {
1402       TrkrDefs::cluskey cluster_key = *iter;
1403       int layer = TrkrDefs::getLayer(cluster_key);
1404       if(layer>6) ntpc++;
1405     }
1406 
1407       if(charge<0) std::cout << "electron: " << pt << " " << eta << " " << nmvtx << " " << ntpc << " " << charge << " " << cemc_ecore << std::endl;
1408       if(charge>0) std::cout << "positron: " << pt << " " << eta << " " << nmvtx << " " << ntpc << " " << charge << " " << cemc_ecore << std::endl;
1409       TLorentzVector tmp = TLorentzVector(px,py,pz,ee);
1410       if(charge>0) {
1411         positrons.push_back(tmp); 
1412         vpeop3x3.push_back(cemc_ecore); vpchi2.push_back(chi2); vpdca.push_back(dca3d_xy);
1413         vpnmvtx.push_back(nmvtx); vpnintt.push_back(nintt); vpntpc.push_back(ntpc);
1414       }
1415       if(charge<0) {
1416         electrons.push_back(tmp); 
1417         vmeop3x3.push_back(cemc_ecore); vmchi2.push_back(chi2); vmdca.push_back(dca3d_xy);
1418         vmnmvtx.push_back(nmvtx); vmnintt.push_back(nintt); vmntpc.push_back(ntpc);
1419       }
1420 
1421   } // end loop over tracks
1422 
1423   double emult = electrons.size();
1424   double pmult = positrons.size();
1425   float tmp1[99];
1426 
1427 // make pairs ----------------------------------------------------------------------------------
1428 
1429   for(long unsigned int i=0; i<electrons.size(); i++) {
1430   for(long unsigned int j=0; j<positrons.size(); j++) {
1431     TLorentzVector pair = electrons[i]+positrons[j];
1432     double mass = pair.M();
1433     hmass->Fill(mass);
1434     cout << "Upsilon mass = " << pair.M() << endl; 
1435     //if(mass<5.0) continue;
1436       tmp1[0] = 1.;
1437       tmp1[1] = pair.M();
1438       tmp1[2] = pair.Pt();
1439       tmp1[3] = pair.Eta();
1440       tmp1[4] = pair.Rapidity();
1441       tmp1[5] = positrons[j].Pt();
1442       tmp1[6] = electrons[i].Pt();
1443       tmp1[7] = positrons[j].Eta();
1444       tmp1[8] = electrons[i].Eta();
1445       tmp1[9] = trackmap->size();
1446       tmp1[10] = emult;
1447       tmp1[11] = pmult;
1448       tmp1[12] = vpchi2[j];
1449       tmp1[13] = vmchi2[i];
1450       tmp1[14] = vpdca[j];
1451       tmp1[15] = vmdca[i];
1452       tmp1[16] = vpeop3x3[j];
1453       tmp1[17] = vmeop3x3[i];
1454       tmp1[18] = electrons[i].DrEtaPhi(positrons[j]);
1455       tmp1[19] = vpnmvtx[j];
1456       tmp1[20] = vmnmvtx[i];
1457       tmp1[21] = vpnintt[j];
1458       tmp1[22] = vmnintt[i];
1459       tmp1[23] = vpntpc[j];
1460       tmp1[24] = vmntpc[i];
1461       tmp1[25] = 0.;
1462        ntp2->Fill(tmp1);
1463   }}
1464 
1465   if(electrons.size()>1) {
1466   for(long unsigned int i=0; i<electrons.size()-1; i++) {
1467   for(long unsigned int j=i+1; j<electrons.size(); j++) {
1468     TLorentzVector pair = electrons[i]+electrons[j];
1469     double mass = pair.M();
1470     hmass->Fill(mass);
1471     //if(mass<5.0) continue;
1472       tmp1[0] = 3.;
1473       tmp1[1] = pair.M();
1474       tmp1[2] = pair.Pt();
1475       tmp1[3] = pair.Eta();
1476       tmp1[4] = pair.Rapidity();
1477       tmp1[5] = electrons[j].Pt();
1478       tmp1[6] = electrons[i].Pt();
1479       tmp1[7] = electrons[j].Eta();
1480       tmp1[8] = electrons[i].Eta();
1481       tmp1[9] = trackmap->size();
1482       tmp1[10] = emult;
1483       tmp1[11] = pmult;
1484       tmp1[12] = vmchi2[j];
1485       tmp1[13] = vmchi2[i];
1486       tmp1[14] = vmdca[j];
1487       tmp1[15] = vmdca[i];
1488       tmp1[16] = vmeop3x3[j];
1489       tmp1[17] = vmeop3x3[i];
1490       tmp1[18] = electrons[i].DrEtaPhi(electrons[j]);
1491       tmp1[19] = vmnmvtx[j];
1492       tmp1[20] = vmnmvtx[i];
1493       tmp1[21] = vmnintt[j];
1494       tmp1[22] = vmnintt[i];
1495       tmp1[23] = vmntpc[j];
1496       tmp1[24] = vmntpc[i];
1497       tmp1[25] = 0.;
1498        ntp2->Fill(tmp1);
1499   }}}
1500 
1501   if(positrons.size()>1) {
1502   for(long unsigned int i=0; i<positrons.size()-1; i++) {
1503   for(long unsigned int j=i+1; j<positrons.size(); j++) {
1504     TLorentzVector pair = positrons[i]+positrons[j];
1505     double mass = pair.M();
1506     hmass->Fill(mass);
1507     //if(mass<5.0) continue;
1508       tmp1[0] = 2.;
1509       tmp1[1] = pair.M();
1510       tmp1[2] = pair.Pt();
1511       tmp1[3] = pair.Eta();
1512       tmp1[4] = pair.Rapidity();
1513       tmp1[5] = positrons[j].Pt();
1514       tmp1[6] = positrons[i].Pt();
1515       tmp1[7] = positrons[j].Eta();
1516       tmp1[8] = positrons[i].Eta();
1517       tmp1[9] = trackmap->size();
1518       tmp1[10] = emult;
1519       tmp1[11] = pmult;
1520       tmp1[12] = vpchi2[j];
1521       tmp1[13] = vpchi2[i];
1522       tmp1[14] = vpdca[j];
1523       tmp1[15] = vpdca[i];
1524       tmp1[16] = vpeop3x3[j];
1525       tmp1[17] = vpeop3x3[i];
1526       tmp1[18] = positrons[i].DrEtaPhi(positrons[j]);
1527       tmp1[19] = vpnmvtx[j];
1528       tmp1[20] = vpnmvtx[i];
1529       tmp1[21] = vpnintt[j];
1530       tmp1[22] = vpnintt[i];
1531       tmp1[23] = vpntpc[j];
1532       tmp1[24] = vpntpc[i];
1533       tmp1[25] = 0.;
1534        ntp2->Fill(tmp1);
1535   }}}
1536 
1537   return Fun4AllReturnCodes::EVENT_OK;
1538 } 
1539 
1540 //======================================================================
1541 
1542 bool sPHAnalysis::isElectron(SvtxTrack* trk)
1543 {
1544   double px = trk->get_px();
1545   double py = trk->get_py();
1546   double pz = trk->get_pz();
1547   double pt = sqrt(px*px+py*py);
1548   double pp = sqrt(pt*pt+pz*pz);
1549   double e3x3 = trk->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::CEMC);
1550   //double eclus = trk->get_cal_cluster_e(SvtxTrack::CAL_LAYER::CEMC);
1551   if(pp==0.) return false;
1552   //if(pt<2.0) return false;
1553   if(pt<0.1) return false;
1554   //cout << e3x3 << " " << eclus << endl;
1555   if(isnan(e3x3)) return false;
1556   //if(e3x3/pp<0.70) return false;
1557   if(e3x3/pp<0.1) return false;
1558   return true;
1559 }
1560 
1561 //======================================================================
1562 
1563 int sPHAnalysis::End(PHCompositeNode *topNode) 
1564 {
1565   std::cout << "sPHAnalysis::End() started, Number of processed events = " << EventNumber << endl;
1566 //   HepMC::write_HepMC_IO_block_end( ascii_io );
1567   std::cout << "Writing out..." << endl;
1568   OutputNtupleFile->Write();
1569   std::cout << "Closing output file..." << endl;
1570   OutputNtupleFile->Close();
1571 //  ofs.close();
1572   return Fun4AllReturnCodes::EVENT_OK;
1573 }
1574 
1575
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!