analysis/PhotonJet/PhotonJet.C

0001 //general fun4all and subsysreco includes
0002 #include <fun4all/Fun4AllServer.h>
0003 #include <g4main/PHG4Hit.h>
0004 #include <g4main/PHG4Particle.h>
0005 #include <g4main/PHG4TruthInfoContainer.h>
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/getClass.h>
0008
0009 //calorimeter includes
0010 #include <calobase/RawCluster.h>
0011 #include <calobase/RawClusterContainer.h>
0012 #include <calobase/RawClusterUtility.h>
0013 #include <calobase/RawTower.h>
0014 #include <calobase/RawTowerContainer.h>
0015 #include <calobase/RawTowerGeom.h>
0016 #include <calobase/RawTowerGeomContainer.h>
0017 #include <calotrigger/CaloTriggerInfo.h>
0018
0019 #include <g4jets/Jet.h>
0020 #include <g4jets/JetMap.h>
0021 #include <g4vertex/GlobalVertex.h>
0022 #include <g4vertex/GlobalVertexMap.h>
0023 #include <trackbase_historic/SvtxTrack.h>
0024 #include <trackbase_historic/SvtxTrackMap.h>
0025 #include <trackbase_historic/SvtxVertex.h>
0026 #include <trackbase_historic/SvtxVertexMap.h>
0027
0028 //evaluation includes
0029 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
0030 #include <g4eval/JetEvalStack.h>
0031 #include <g4eval/SvtxEvalStack.h>
0032
0033 //hepmc includes
0034 #include <HepMC/GenEvent.h>
0035 #include <HepMC/GenVertex.h>
0036 #include <phhepmc/PHHepMCGenEvent.h>
0037 #include <phhepmc/PHHepMCGenEventMap.h>
0038
0039 //general c++ includes
0040 #include <TLorentzVector.h>
0041 #include <cassert>
0042 #include <iostream>
0043 #include <sstream>
0044
0045 #include "PhotonJet.h"
0046
0047 using namespace std;
0048
0049 PhotonJet::PhotonJet(const std::string &name)
0050   : SubsysReco("PHOTONJET")
0051 {
0052   //just set all global values to -999 from the start so that they have a spot in memory
0053   initialize_values();
0054
0055   outfilename = name;
0056
0057   //initialize public member values
0058
0059   //default don't use isocone algorithm
0060   use_isocone = 0;
0061
0062   //default do not use tracked jets (or tracks)
0063   eval_tracked_jets = 0;
0064
0065   //default use 0.3 jet cone
0066   jet_cone_size = 3;
0067
0068   //default set beginning number of events to 0
0069   //this can be changed with one of the public functions if e.g. you want individual event number IDs for multiple MC blocks of events
0070   nevents = 0;
0071
0072   //default to barrel sPHENIX acceptance
0073   etalow = -1;
0074   etahigh = 1;
0075
0076   //default jetminptcut
0077   minjetpt = 5.;
0078
0079   //default mincluscut
0080   mincluspt = 0.5;
0081
0082   //default minimum direct photon p_T
0083   mindp_pt = 10;
0084
0085   //default to use the trigger emulator
0086   usetrigger = 1;
0087
0088   //default to using position corrected emcal clusters
0089   use_pos_cor_cemc = 1;
0090
0091   //default to not AA and not using embedded background subtraction
0092   is_AA = 0;
0093 }
0094
0095 int PhotonJet::Init(PHCompositeNode *topnode)
0096 {
0097   if (Verbosity() > 1)
0098   {
0099     cout << "COLLECTING PHOTON-JET PAIRS FOR THE FOLLOWING: " << endl;
0100     cout << "GATHERING JETS: " << jet_cone_size << endl;
0101     cout << "GATHERING IN ETA: [" << etalow
0102          << "," << etahigh << "]" << endl;
0103     cout << "EVALUATING TRACKED JETS: " << eval_tracked_jets << endl;
0104     cout << "USING ISOLATION CONE: " << use_isocone << endl;
0105   }
0106
0107   //create output photonjet tfile which contains output TTrees
0108   file = new TFile(outfilename.c_str(), "RECREATE");
0109
0110   //create some basic histograms
0111   ntruthconstituents_h = new TH1F("ntruthconstituents", "", 200, 0, 200);
0112   percent_photon_h = new TH1F("percent_photon_h",
0113                               ";E_{photon}/E_{jet}; Counts", 200, 0, 2);
0114   histo = new TH1F("histo", "histo", 100, -3, 3);
0115
0116   //create basic tree
0117   tree = new TTree("tree", "a tree");
0118   tree->Branch("nevents", &nevents, "nevents/I");
0119
0120   //main trees are set in this fxn - branch addresses are defined to global data types
0121   Set_Tree_Branches();
0122
0123   return 0;
0124 }
0125
0126 int PhotonJet::process_event(PHCompositeNode *topnode)
0127 {
0128   // event number tracker,
0129   if (nevents % 10 == 0)
0130     cout << "at event number " << nevents << endl;
0131
0132   //get the requested size jets
0133   ostringstream truthjetsize;
0134   ostringstream recojetsize;
0135   ostringstream trackjetsize;
0136
0137   //these are the node names for Truth, Calorimeter tower, and tracked jets
0138   truthjetsize.str("");
0139   truthjetsize << "AntiKt_Truth_r";
0140   recojetsize.str("");
0141   recojetsize << "AntiKt_Tower_r";
0142   trackjetsize.str("");
0143   trackjetsize << "AntiKt_Track_r";
0144
0145   if (jet_cone_size == 2)
0146   {
0147     truthjetsize << "02";
0148     recojetsize << "02";
0149     trackjetsize << "02";
0150   }
0151   else if (jet_cone_size == 3)
0152   {
0153     truthjetsize << "03";
0154     recojetsize << "03";
0155     trackjetsize << "03";
0156   }
0157   else if (jet_cone_size == 4)
0158   {
0159     truthjetsize << "04";
0160     recojetsize << "04";
0161     trackjetsize << "04";
0162   }
0163   else if (jet_cone_size == 5)
0164   {
0165     truthjetsize << "05";
0166     recojetsize << "05";
0167     trackjetsize << "05";
0168   }
0169   else if (jet_cone_size == 6)
0170   {
0171     truthjetsize << "06";
0172     recojetsize << "06";
0173     trackjetsize << "06";
0174   }
0175
0176   else if (jet_cone_size == 7)
0177   {
0178     truthjetsize << "07";
0179     recojetsize << "07";
0180     trackjetsize << "07";
0181   }
0182   //if its some other number just set it to 0.4
0183
0184   else
0185   {
0186     truthjetsize << "04";
0187     recojetsize << "04";
0188   }
0189
0190   if (Verbosity() > 1)
0191   {
0192     cout << "Gathering RECO Jets:  " << recojetsize.str().c_str() << endl;
0193     cout << "Gathering TRUTH Jets:  " << truthjetsize.str().c_str() << endl;
0194   }
0195
0196   //get the nodes from the NodeTree
0197
0198   //JetMap nodes
0199   JetMap *truth_jets =
0200       findNode::getClass<JetMap>(topnode, truthjetsize.str().c_str());
0201
0202   JetMap *reco_jets = 0;
0203   if (!is_AA)
0204   {
0205     reco_jets = findNode::getClass<JetMap>(topnode, recojetsize.str().c_str());
0206   }
0207
0208   JetMap *tracked_jets =
0209       findNode::getClass<JetMap>(topnode, trackjetsize.str().c_str());
0210
0211   //G4 truth particle node
0212   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topnode, "G4TruthInfo");
0213
0214   //EMCal raw tower cluster node
0215   RawClusterContainer *clusters = 0;
0216   if (!use_pos_cor_cemc)
0217     clusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_CEMC");
0218
0219   //EMCal position calibrated cluster node
0220   if (use_pos_cor_cemc)
0221     clusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_POS_COR_CEMC");
0222
0223   //SVTX tracks node
0224   SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topnode, "SvtxTrackMap");
0225
0226   //trigger emulator
0227   CaloTriggerInfo *trigger = findNode::getClass<CaloTriggerInfo>(topnode, "CaloTriggerInfo");
0228
0229   PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0230
0231   //for truth track matching
0232   SvtxEvalStack *svtxevalstack = new SvtxEvalStack(topnode);
0233   svtxevalstack->next_event(topnode);
0234
0235   if (is_AA)
0236   {
0237     recojetsize << "_Sub1";
0238     reco_jets =
0239         findNode::getClass<JetMap>(topnode, recojetsize.str().c_str());
0240   }
0241
0242   //for truth jet matching
0243   JetEvalStack *_jetevalstack = 0;
0244   if (!is_AA)
0245     _jetevalstack =
0246         new JetEvalStack(topnode, recojetsize.str().c_str(),
0247                          truthjetsize.str().c_str());
0248   //the jet eval stack doesn't work for AA - so this is useless
0249   //in that case the code is setup to match reco and truth jets based on
0250   //their difference in deltaphi and deltaeta
0251   else
0252     _jetevalstack = new JetEvalStack(topnode,
0253                                      "AntiKt_Tower_r04_Sub1",
0254                                      "AntiKt_Truth_r04");
0255
0256   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topnode, "GlobalVertexMap");
0257   if (!vertexmap)
0258   {
0259     cout << "PhotonJet::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << endl;
0260     assert(vertexmap);  // force quit
0261
0262     return 0;
0263   }
0264
0265   if (vertexmap->empty())
0266   {
0267     cout << "PhotonJet::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << endl;
0268     return 0;
0269   }
0270
0271   GlobalVertex *vtx = vertexmap->begin()->second;
0272   if (vtx == nullptr) return 0;
0273
0274
0275   RawTowerContainer *_cemctowers = findNode::getClass<RawTowerContainer>(topnode,"TOWER_CALIB_CEMC");
0276   RawTowerContainer *_hcalintowers = findNode::getClass<RawTowerContainer>(topnode,"TOWER_CALIB_HCALIN");
0277   RawTowerContainer *_hcalouttowers = findNode::getClass<RawTowerContainer>(topnode,"TOWER_CALIB_HCALOUT");
0278   RawTowerGeomContainer *_cemctowergeom = findNode::getClass<RawTowerGeomContainer>(topnode,"TOWERGEOM_CEMC");
0279   RawTowerGeomContainer *_hcalintowergeom = findNode::getClass<RawTowerGeomContainer>(topnode,"TOWERGEOM_HCALIN");
0280   RawTowerGeomContainer *_hcalouttowergeom = findNode::getClass<RawTowerGeomContainer>(topnode,"TOWERGEOM_HCALOUT");
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290   //Make sure all nodes for analysis are here. If one isn't in the NodeTree, bail
0291   if (!trigger && usetrigger)
0292   {
0293     cout << "NO TRIGGER EMULATOR, BAILING" << endl;
0294     return 0;
0295   }
0296   if (!tracked_jets && eval_tracked_jets)
0297   {
0298     cout << "no tracked jets, bailing" << endl;
0299     return 0;
0300   }
0301   if (!truth_jets)
0302   {
0303     cout << "no truth jets" << endl;
0304     return 0;
0305   }
0306   if (!reco_jets)
0307   {
0308     cout << "no reco jets" << endl;
0309     return 0;
0310   }
0311   if (!truthinfo)
0312   {
0313     cout << "no truth track info" << endl;
0314     return 0;
0315   }
0316   if (!clusters)
0317   {
0318     cout << "no cluster info" << endl;
0319     return 0;
0320   }
0321   if (!trackmap && eval_tracked_jets)
0322   {
0323     cout << "no track map" << endl;
0324     return 0;
0325   }
0326   if (!_jetevalstack)
0327   {
0328     cout << "no good truth jets" << endl;
0329     return 0;
0330   }
0331
0332   //For jet and track truth/reco matching
0333   JetRecoEval *recoeval = _jetevalstack->get_reco_eval();
0334   SvtxTrackEval *trackeval = svtxevalstack->get_track_eval();
0335   JetTruthEval *trutheval = _jetevalstack->get_truth_eval();
0336
0337
0338
0339   //now we have all the nodes, so collect the data from the various nodes
0340
0341   /***********************************************
0342
0343   GET ALL THE HEPMC EVENT LEVEL TRUTH PARTICLES
0344
0345   ************************************************/
0346
0347   PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topnode, "PHHepMCGenEventMap");
0348
0349   if (!hepmceventmap)
0350   {
0351     cout << "hepmc event map node is missing, BAILING" << endl;
0352     //return 0;
0353   }
0354
0355   if (Verbosity() > 1)
0356   {
0357     cout << "Getting HEPMC truth particles " << endl;
0358   }
0359
0360   //you can iterate over the number of events in a hepmc event
0361   //for pile up events where you have multiple hard scatterings per bunch crossing
0362   for (PHHepMCGenEventMap::ConstIter iterr = hepmceventmap->begin();
0363        iterr != hepmceventmap->end();
0364        ++iterr)
0365   {
0366     PHHepMCGenEvent *hepmcevent = iterr->second;
0367
0368     if (hepmcevent)
0369     {
0370       //get the event characteristics, inherited from HepMC classes
0371       HepMC::GenEvent *truthevent = hepmcevent->getEvent();
0372       if (!truthevent)
0373       {
0374         cout << PHWHERE << "no evt pointer under phhepmvgeneventmap found " << endl;
0375         return 0;
0376       }
0377
0378       //get the interacting protons
0379       if (truthevent->valid_beam_particles())
0380       {
0381         HepMC::GenParticle *part1 = truthevent->beam_particles().first;
0382
0383         //beam energy
0384         beam_energy = fabs(part1->momentum().pz());
0385       }
0386
0387       //Parton info
0388       HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
0389
0390       partid1 = pdfinfo->id1();
0391       partid2 = pdfinfo->id2();
0392       x1 = pdfinfo->x1();
0393       x2 = pdfinfo->x2();
0394
0395       //are there multiple partonic intercations in a p+p event
0396       mpi = truthevent->mpi();
0397
0398       numparticlesinevent = 0;
0399
0400       //Get the PYTHIA signal process id identifying the 2-to-2 hard process
0401       process_id = truthevent->signal_process_id();
0402
0403       //loop over all the truth particles and get their information
0404       for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
0405            iter != truthevent->particles_end();
0406            ++iter)
0407       {
0408         //get each pythia particle characteristics
0409         truthenergy = (*iter)->momentum().e();
0410         truthpid = (*iter)->pdg_id();
0411
0412         trutheta = (*iter)->momentum().pseudoRapidity();
0413         truthphi = (*iter)->momentum().phi();
0414         truthpx = (*iter)->momentum().px();
0415         truthpy = (*iter)->momentum().py();
0416         truthpz = (*iter)->momentum().pz();
0417         truthpt = sqrt(truthpx * truthpx + truthpy * truthpy);
0418
0419         //fill the truth tree
0420         truthtree->Fill();
0421         numparticlesinevent++;
0422       }
0423     }
0424   }
0425
0426   //these are global variables, i.e. the highest pt cluster in an event
0427   cluseventenergy = 0;
0428   cluseventphi = 0;
0429   cluseventpt = 0;
0430   cluseventeta = 0;
0431   float lastenergy = 0;
0432
0433   if (Verbosity() > 1)
0434   {
0435     cout << "get G4 stable truth particles" << endl;
0436   }
0437
0438   //loop over the G4 truth (stable) particles
0439   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0440        iter != range.second;
0441        ++iter)
0442   {
0443     //get this particle
0444     PHG4Particle *truth = iter->second;
0445
0446     //get this particles momentum, etc.
0447     truthpx = truth->get_px();
0448     truthpy = truth->get_py();
0449     truthpz = truth->get_pz();
0450     truthp = sqrt(truthpx * truthpx + truthpy * truthpy + truthpz * truthpz);
0451     truthenergy = truth->get_e();
0452
0453     TLorentzVector vec;
0454     vec.SetPxPyPzE(truthpx, truthpy, truthpz, truthenergy);
0455
0456     truthpt = sqrt(truthpx * truthpx + truthpy * truthpy);
0457
0458     truthphi = vec.Phi();
0459
0460     trutheta = TMath::ATanH(truthpz / truthenergy);
0461     //check for nans
0462     if (trutheta != trutheta)
0463       trutheta = -9999;
0464     truthpid = truth->get_pid();
0465
0466     //find the highest energy photon in the event in the eta range
0467     if (truthpid == 22 && truthenergy > lastenergy
0468     && trutheta > etalow && trutheta < etahigh)
0469     {
0470       lastenergy = truthenergy;
0471       cluseventenergy = truthenergy;
0472       cluseventpt = truthpt;
0473       cluseventphi = truthphi;
0474       cluseventeta = trutheta;
0475     }
0476     //fill the g4 truth tree
0477     truth_g4particles->Fill();
0478   }
0479
0480   /***************************************
0481
0482    TRUTH JETS
0483
0484    ***************************************/
0485   if (Verbosity() > 1)
0486   {
0487     cout << "get the truth jets" << endl;
0488   }
0489
0490   //loop over the truth jets
0491   float hardest_jet = 0;
0492   // If not truth jet is found that satisfies the requirements,
0493   // then the jet 4 vector will be (0,0,0,0)
0494   hardest_jetpt = 0;
0495   hardest_jetenergy = 0;
0496   hardest_jeteta = 0;
0497   hardest_jetphi = 0;
0498   for (JetMap::Iter iter = truth_jets->begin();
0499        iter != truth_jets->end();
0500        ++iter)
0501   {
0502     Jet *jet = iter->second;
0503
0504     truthjetpt = jet->get_pt();
0505
0506     //only collect truthjets above the minjetpt cut
0507     if (truthjetpt < minjetpt)
0508       continue;
0509
0510     truthjeteta = jet->get_eta();
0511
0512     //make the width extra just to be safe and collect truth jets
0513     //that might be e.g. half in the sphenix acceptance
0514     if (truthjeteta < (etalow - 1) || truthjeteta > (etahigh + 1))
0515       continue;
0516
0517     truthjetpx = jet->get_px();
0518     truthjetpy = jet->get_py();
0519     truthjetpz = jet->get_pz();
0520     truthjetphi = jet->get_phi();
0521     truthjetmass = jet->get_mass();
0522     truthjetp = jet->get_p();
0523     truthjetenergy = jet->get_e();
0524
0525     //check that the jet isn't just a high pt photon
0526
0527     //get the truth constituents of the jet
0528     std::set<PHG4Particle *> truthjetcomp =
0529         trutheval->all_truth_particles(jet);
0530     ntruthconstituents = 0;
0531     float truthjetenergysum = 0;
0532     float truthjethighestphoton = 0;
0533     //loop over the constituents of the truth jet
0534     for (std::set<PHG4Particle *>::iterator iter2 = truthjetcomp.begin();
0535          iter2 != truthjetcomp.end();
0536          ++iter2)
0537     {
0538       //get the particle of the truthjet
0539       PHG4Particle *truthpart = *iter2;
0540       if (!truthpart)
0541       {
0542         cout << "no truth particles in the jet??" << endl;
0543         break;
0544       }
0545
0546       ntruthconstituents++;
0547
0548       int constpid = truthpart->get_pid();
0549       float conste = truthpart->get_e();
0550
0551       truthjetenergysum += conste;
0552
0553       if (constpid == 22)
0554       {
0555         if (conste > truthjethighestphoton)
0556           truthjethighestphoton = conste;
0557       }
0558     }
0559     ntruthconstituents_h->Fill(ntruthconstituents);
0560
0561     //we want jets that are real jets, not just an e.g. isolated photon
0562     //require that the number of constituents in the jet be larger than 3
0563     float percent_photon = truthjethighestphoton / truthjetenergy;
0564     percent_photon_h->Fill(percent_photon);
0565
0566     //if there is a high energy photon that is 80% of the jets energy, skip it
0567     //it is likely the near-side direct photon
0568     if (percent_photon > 0.8)
0569     {
0570       continue;
0571     }
0572
0573     //we also want jets to have at least 3 constituents
0574     if (ntruthconstituents < 3)
0575       continue;
0576
0577     //identify the highest pt jet in the event
0578     if (truthjetpt > hardest_jet)
0579     {
0580       hardest_jet = truthjetpt;
0581       hardest_jetpt = truthjetpt;
0582       hardest_jetenergy = truthjetenergy;
0583       hardest_jeteta = truthjeteta;
0584       hardest_jetphi = truthjetphi;
0585     }
0586
0587     //fill the truthjet tree
0588     truthjettree->Fill();
0589   }
0590
0591   //fill the event tree with e.g. x1,x2 partonic momentum fractions,
0592   //process id, highest energy photon, etc.
0593   event_tree->Fill();
0594
0595   if (Verbosity() > 1)
0596   {
0597     cout << "get trigger emulator info" << endl;
0598   }
0599   /***********************************************
0600
0601    TRIGGER EMULATOR INFO
0602
0603   ************************************************/
0604   if (usetrigger)
0605   {
0606     //get the 4x4 trigger tile info
0607     E_4x4 = trigger->get_best_EMCal_4x4_E();
0608     phi_4x4 = trigger->get_best_EMCal_4x4_phi();
0609     eta_4x4 = trigger->get_best_EMCal_4x4_eta();
0610
0611     //get the 2x2 trigger tile info
0612     E_2x2 = trigger->get_best_EMCal_2x2_E();
0613     phi_2x2 = trigger->get_best_EMCal_2x2_phi();
0614     eta_2x2 = trigger->get_best_EMCal_2x2_eta();
0615   }
0616
0617   /***********************************************
0618
0619   GET THE EMCAL CLUSTERS
0620
0621   ************************************************/
0622   if (Verbosity() > 1)
0623   {
0624     cout << "Get EMCal Cluster" << endl;
0625   }
0626
0627   RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0628   RawClusterContainer::ConstIterator clusiter;
0629
0630   //loop over the emcal clusters
0631   for (clusiter = begin_end.first;
0632        clusiter != begin_end.second;
0633        ++clusiter)
0634   {
0635     //get this cluster
0636     RawCluster *cluster = clusiter->second;
0637
0638     //get cluster characteristics
0639     //this helper class determines the photon characteristics
0640     //depending on the vertex position
0641     //this is important for e.g. eta determination and E_T determination
0642     CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0643     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
0644     clus_energy = E_vec_cluster.mag();
0645     clus_eta = E_vec_cluster.pseudoRapidity();
0646     clus_theta = E_vec_cluster.getTheta();
0647     clus_pt = E_vec_cluster.perp();
0648     clus_phi = E_vec_cluster.getPhi();
0649
0650     if (clus_pt < mincluspt)
0651       continue;
0652
0653     if (clus_eta < etalow)
0654       continue;
0655     if (clus_eta > etahigh)
0656       continue;
0657
0658     TLorentzVector *clus = new TLorentzVector();
0659     clus->SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clus_energy);
0660
0661     float dumarray[4];
0662     clus->GetXYZT(dumarray);
0663     clus_x = dumarray[0];
0664     clus_y = dumarray[1];
0665     clus_z = dumarray[2];
0666     clus_t = dumarray[3];
0667
0668     clus_px = clus_pt * TMath::Cos(clus_phi);
0669     clus_py = clus_pt * TMath::Sin(clus_phi);
0670     clus_pz = sqrt(clus_energy * clus_energy - clus_px * clus_px - clus_py * clus_py);
0671
0672     //fill the cluster tree with all emcal clusters
0673     cluster_tree->Fill();
0674
0675     //now determine the direct photon
0676     //only interested in high pt photons to be isolated i.e. direct photons
0677     if (clus_pt < mindp_pt)
0678       continue;
0679     //require that the entire isolation cone fall within sPHENIX acceptance
0680     if (fabs(clus_eta) > (1.0 - isoconeradius) && use_isocone)
0681       continue;
0682
0683     if (use_isocone)
0684     {
0685       //get the energy sum in the cone surrounding the photon
0686       float energysum = ConeSum(cluster, clusters, trackmap, isoconeradius, vtx);
0687
0688       //check if energy sum is less than 10% of isolated photon energy
0689       bool conecut = energysum > 0.1 * clus_energy;
0690       if (conecut)
0691         continue;
0692     }
0693
0694     //find the associated truth high pT photon with this reconstructed photon
0695     for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0696          iter != range.second;
0697          ++iter)
0698     {
0699       //get this truth particle
0700       PHG4Particle *truth = iter->second;
0701
0702       clustruthpid = truth->get_pid();
0703
0704       //check that it is a photon, i.e. pid==22
0705       if (clustruthpid == 22)
0706       {
0707         clustruthpx = truth->get_px();
0708         clustruthpy = truth->get_py();
0709         clustruthpz = truth->get_pz();
0710         clustruthenergy = truth->get_e();
0711         clustruthpt = sqrt(clustruthpx * clustruthpx + clustruthpy * clustruthpy);
0712         if (clustruthpt < mindp_pt)
0713           continue;
0714
0715         TLorentzVector vec;
0716         vec.SetPxPyPzE(clustruthpx, clustruthpy, clustruthpz, clustruthenergy);
0717         clustruthphi = vec.Phi();
0718
0719         clustrutheta = TMath::ATanH(clustruthpz / clustruthenergy);
0720         //check for nans
0721         if (clustrutheta != clustrutheta)
0722           clustrutheta = -9999;
0723
0724         //check that the truth photon has similar eta/phi to reco photon
0725         //the clustering has a resolution of about 0.005 rads so this is sufficient
0726         if (fabs(clustruthphi - clus_phi) > 0.02 || fabs(clustrutheta - clus_eta) > 0.02)
0727           continue;
0728
0729         //once the photon is found and the values are set
0730         //just break out of the loop
0731         break;
0732       }
0733     }
0734     //fill isolated clusters tree, i.e. all isolated clusters regardless
0735     //of if an away-side jet is found also
0736     isolated_clusters->Fill();
0737
0738     //get back to back reconstructed hadrons/jets for photon-jet processes
0739     //two different functions for Au+Au vs p+p due to the way truth jet matching
0740     //is performed between the two systems
0741     if (!is_AA)
0742       GetRecoHadronsAndJets(cluster, trackmap,
0743                             reco_jets, tracked_jets,
0744                             recoeval, trackeval,
0745                             truthinfo,
0746                             trutheval,
0747                             truth_jets,
0748                             vtx);
0749
0750     if (is_AA)
0751       GetRecoHadronsAndJetsAA(cluster,
0752                               trackmap,
0753                               reco_jets,
0754                               truthinfo,
0755                               truth_jets,
0756                               vtx);
0757   }
0758
0759   /***********************************************
0760
0761   GET THE SVTX TRACKS
0762
0763   ************************************************/
0764   if (eval_tracked_jets)
0765   {
0766     if (Verbosity() > 1)
0767     {
0768       cout << "Get the Tracks" << endl;
0769     }
0770     for (SvtxTrackMap::Iter iter = trackmap->begin();
0771          iter != trackmap->end();
0772          ++iter)
0773     {
0774       SvtxTrack *track = iter->second;
0775
0776       //get reco info
0777       tr_px = track->get_px();
0778       tr_py = track->get_py();
0779       tr_pz = track->get_pz();
0780       tr_p = sqrt(tr_px * tr_px + tr_py * tr_py + tr_pz * tr_pz);
0781
0782       tr_pt = sqrt(tr_px * tr_px + tr_py * tr_py);
0783
0784       if (tr_pt < 0.5)
0785         continue;
0786
0787       tr_phi = track->get_phi();
0788       tr_eta = track->get_eta();
0789
0790       if (tr_eta < etalow || tr_eta > etahigh)
0791         continue;
0792
0793       charge = track->get_charge();
0794       chisq = track->get_chisq();
0795       ndf = track->get_ndf();
0796       dca = track->get_dca();
0797       tr_x = track->get_x();
0798       tr_y = track->get_y();
0799       tr_z = track->get_z();
0800
0801       //get truth track info
0802       PHG4Particle *truthtrack = trackeval->max_truth_particle_by_nclusters(track);
0803       truth_is_primary = truthinfo->is_primary(truthtrack);
0804
0805       truthtrackpx = truthtrack->get_px();
0806       truthtrackpy = truthtrack->get_py();
0807       truthtrackpz = truthtrack->get_pz();
0808       truthtrackp = sqrt(truthtrackpx * truthtrackpx + truthtrackpy * truthtrackpy + truthtrackpz * truthtrackpz);
0809
0810       truthtracke = truthtrack->get_e();
0811
0812       TLorentzVector *Truthtrack = new TLorentzVector();
0813       Truthtrack->SetPxPyPzE(truthtrackpx, truthtrackpy, truthtrackpz, truthtracke);
0814
0815       truthtrackpt = Truthtrack->Pt();
0816       truthtrackphi = Truthtrack->Phi();
0817       truthtracketa = Truthtrack->Eta();
0818       truthtrackpid = truthtrack->get_pid();
0819
0820       tracktree->Fill();
0821     }
0822   }
0823
0824   /***************************************
0825
0826    RECONSTRUCTED JETS
0827
0828    ***************************************/
0829   if (Verbosity() > 1)
0830   {
0831     cout << "Get all Reco Jets" << endl;
0832   }
0833
0834   for (JetMap::Iter iter = reco_jets->begin();
0835        iter != reco_jets->end();
0836        ++iter)
0837   {
0838     Jet *jet = iter->second;
0839     Jet *truthjet = recoeval->max_truth_jet_by_energy(jet);
0840     recojetpt = jet->get_pt();
0841     if (recojetpt < minjetpt)
0842       continue;
0843
0844     recojeteta = jet->get_eta();
0845
0846     //reco jet eta better not be outside of the sPHENIX acceptance....
0847     if (recojeteta < (etalow - 1) || recojeteta > (etahigh + 1))
0848       continue;
0849
0850     //get reco jet characteristics
0851     recojetid = jet->get_id();
0852     recojetpx = jet->get_px();
0853     recojetpy = jet->get_py();
0854     recojetpz = jet->get_pz();
0855     recojetphi = jet->get_phi();
0856     recojetmass = jet->get_mass();
0857     recojetp = jet->get_p();
0858     recojetenergy = jet->get_e();
0859
0860     //if truthjet exists, then it is p+p and we can use the stackeval
0861     //for truthjet matching
0862     if (truthjet)
0863     {
0864       truthjetid = truthjet->get_id();
0865       truthjetp = truthjet->get_p();
0866       truthjetphi = truthjet->get_phi();
0867       truthjeteta = truthjet->get_eta();
0868       truthjetpt = truthjet->get_pt();
0869       truthjetenergy = truthjet->get_e();
0870       truthjetpx = truthjet->get_px();
0871       truthjetpy = truthjet->get_py();
0872       truthjetpz = truthjet->get_pz();
0873       truthjetmass = truthjet->get_mass();
0874
0875       //check that the jet isn't just a photon or something like this
0876       //needs at least 2 constituents and that 80% of jet isn't from one photon
0877       std::set<PHG4Particle *> truthjetcomp =
0878           trutheval->all_truth_particles(truthjet);
0879       ntruthconstituents = 0;
0880       float truthjetenergysum = 0;
0881       float truthjethighestphoton = 0;
0882       for (std::set<PHG4Particle *>::iterator iter2 = truthjetcomp.begin();
0883            iter2 != truthjetcomp.end();
0884            ++iter2)
0885       {
0886         PHG4Particle *truthpart = *iter2;
0887         if (!truthpart)
0888         {
0889           cout << "no truth particles in the jet??" << endl;
0890           break;
0891         }
0892         ntruthconstituents++;
0893
0894         int constpid = truthpart->get_pid();
0895         float conste = truthpart->get_e();
0896
0897         truthjetenergysum += conste;
0898
0899         if (constpid == 22)
0900         {
0901           if (conste > truthjethighestphoton)
0902             truthjethighestphoton = conste;
0903     }
0904       }
0905
0906       //if the highest energy photon in the jet is 80% of the jets energy
0907       //its probably an isolated photon and so we want to not include it in the tree
0908       float percent_photon = truthjethighestphoton / truthjetenergy;
0909       if (percent_photon > 0.8)
0910       {
0911         continue;
0912       }
0913       if (!is_AA && ntruthconstituents < 3)
0914         continue;
0915     }
0916
0917     //if truthjet was null then we just match the reco-truth jets by
0918     //their distance in dphi deta space
0919     else
0920     {
0921       if (Verbosity() > 1)
0922       {
0923         cout << "matching by distance jet" << endl;
0924       }
0925
0926       truthjetid = 0;
0927       truthjetp = 0;
0928       truthjetphi = 0;
0929       truthjeteta = 0;
0930       truthjetpt = 0;
0931       truthjetenergy = 0;
0932       truthjetpx = 0;
0933       truthjetpy = 0;
0934       truthjetpz = 0;
0935
0936       //if the jetevalstack can't find a good truth jet
0937       //try to match reco jet with closest truth jet
0938       float closestjet = 9999;
0939       for (JetMap::Iter iter3 = truth_jets->begin();
0940            iter3 != truth_jets->end();
0941            ++iter3)
0942       {
0943         Jet *jet = iter3->second;
0944
0945         float thisjetpt = jet->get_pt();
0946         if (thisjetpt < minjetpt)
0947           continue;
0948
0949         float thisjeteta = jet->get_eta();
0950         float thisjetphi = jet->get_phi();
0951
0952         float dphi = recojetphi - thisjetphi;
0953         if (dphi > 3. * pi / 2.)
0954           dphi -= pi * 2.;
0955         if (dphi < -1. * pi / 2.)
0956           dphi += pi * 2.;
0957
0958         float deta = recojeteta - thisjeteta;
0959
0960         dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
0961
0962         if (dR < reco_jets->get_par() && dR < closestjet)
0963         {
0964           truthjetid = -9999;  //indicates matched with distance, not jetevalstack
0965           truthjetp = jet->get_p();
0966           truthjetphi = jet->get_phi();
0967           truthjeteta = jet->get_eta();
0968           truthjetpt = jet->get_pt();
0969           truthjetenergy = jet->get_e();
0970           truthjetpx = jet->get_px();
0971           truthjetpy = jet->get_py();
0972           truthjetpz = jet->get_pz();
0973           truthjetmass = jet->get_mass();
0974         }
0975
0976         if (dR < closestjet)
0977         {
0978           closestjet = dR;
0979         }
0980       }
0981     }
0982
0983     //get the reco jet constituents and calculate the dphi deta
0984     //from the jet axis. Added to understand the truth jet - reco jet
0985     //azimuthal offset
0986
0987
0988
0989     for(Jet::ConstIter constiter = jet->begin_comp();
0990     constiter != jet->end_comp();
0991     ++constiter)
0992       {
0993     Jet::SRC source = constiter->first;
0994     unsigned int index = constiter->second;
0995
0996     RawTower *thetower = 0;
0997     if(source == Jet::CEMC_TOWER)
0998       {
0999         thetower = _cemctowers->getTower(index);
1000       }
1001     else if(source == Jet::HCALIN_TOWER)
1002       {
1003         thetower = _hcalintowers->getTower(index);
1004       }
1005     else if(source == Jet::HCALOUT_TOWER)
1006       {
1007         thetower = _hcalouttowers->getTower(index);
1008       }
1009
1010     assert(thetower);
1011
1012     int tower_phi_bin = thetower->get_binphi();
1013     int tower_eta_bin = thetower->get_bineta();
1014
1015     double constphi = -9999;
1016     double consteta = -9999;
1017     if(source == Jet::CEMC_TOWER)
1018       {
1019         constphi = _cemctowergeom->get_phicenter(tower_phi_bin);
1020         consteta = _cemctowergeom->get_etacenter(tower_eta_bin);
1021       }
1022     else if(source == Jet::HCALIN_TOWER)
1023       {
1024         constphi = _hcalintowergeom->get_phicenter(tower_phi_bin);
1025         consteta = _hcalintowergeom->get_etacenter(tower_eta_bin);
1026       }
1027     else if(source == Jet::HCALOUT_TOWER)
1028       {
1029         constphi = _hcalouttowergeom->get_phicenter(tower_phi_bin);
1030         consteta = _hcalouttowergeom->get_etacenter(tower_eta_bin);
1031       }
1032     float checkdphi = truthjetphi - constphi;
1033     if(checkdphi < -1 * TMath::Pi() / 2.)
1034       checkdphi += 2. * TMath::Pi();
1035     if(checkdphi > 3. * TMath::Pi() / 2.)
1036       checkdphi -= 2. * TMath::Pi();
1037
1038     constituent_dphis.push_back(checkdphi);
1039     constituent_detas.push_back(truthjeteta - consteta);
1040
1041       }
1042
1043     recojettree->Fill();
1044
1045     constituent_dphis.resize(0);
1046     constituent_detas.resize(0);
1047
1048   }
1049
1050   if (Verbosity() > 1)
1051   {
1052     cout << "finished event" << endl;
1053   }
1054
1055   nevents++;
1056   tree->Fill();
1057   return 0;
1058 }
1059
1060 int PhotonJet::End(PHCompositeNode *topnode)
1061 {
1062   std::cout << " DONE PROCESSING PHOTONJET PACKAGE" << endl;
1063
1064   file->Write();
1065   file->Close();
1066   return 0;
1067 }
1068 void PhotonJet::GetRecoHadronsAndJetsAA(RawCluster *trig,
1069                                         SvtxTrackMap *tracks,
1070                                         JetMap *recojets,
1071                                         PHG4TruthInfoContainer *alltruth,
1072                                         JetMap *truthjets,
1073                                         GlobalVertex *vtx)
1074 {
1075   CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
1076   CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*trig, vertex);
1077
1078   float trig_phi = E_vec_cluster.getPhi();
1079   float trig_eta = E_vec_cluster.pseudoRapidity();
1080
1081   PHG4TruthInfoContainer::Range range = alltruth->GetPrimaryParticleRange();
1082
1083   //find the matching truth photon to the reco photon
1084   //for the values in the photon-jet tree
1085   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
1086        iter != range.second;
1087        ++iter)
1088   {
1089     PHG4Particle *truth = iter->second;
1090
1091     clustruthpid = truth->get_pid();
1092     if (clustruthpid == 22)
1093     {
1094       clustruthpx = truth->get_px();
1095       clustruthpy = truth->get_py();
1096       clustruthpz = truth->get_pz();
1097       clustruthenergy = truth->get_e();
1098       clustruthpt = sqrt(clustruthpx * clustruthpx + clustruthpy * clustruthpy);
1099       if (clustruthpt < mincluspt)
1100         continue;
1101
1102       TLorentzVector vec;
1103       vec.SetPxPyPzE(clustruthpx, clustruthpy, clustruthpz, clustruthenergy);
1104       clustruthphi = vec.Phi();
1105
1106       clustrutheta = TMath::ATanH(clustruthpz / clustruthenergy);
1107       //check for nans
1108       if (clustrutheta != clustrutheta)
1109         clustrutheta = -9999;
1110
1111       if (fabs(clustruthphi - trig_phi) > 0.02 || fabs(clustrutheta - trig_eta) > 0.02)
1112         continue;
1113
1114       //once the values are set, we've found the truth photon so just break out of this loop
1115       break;
1116     }
1117   }
1118
1119   //find the away-side jets from the direct photon
1120   for (JetMap::Iter iter = recojets->begin();
1121        iter != recojets->end();
1122        ++iter)
1123   {
1124     Jet *jet = iter->second;
1125
1126     _recojetpt = jet->get_pt();
1127     if (_recojetpt < minjetpt)
1128       continue;
1129
1130     _recojeteta = jet->get_eta();
1131     if (_recojeteta < etalow || _recojeteta > etahigh)
1132       continue;
1133
1134     _recojetpx = jet->get_px();
1135     _recojetpy = jet->get_py();
1136     _recojetpz = jet->get_pz();
1137     _recojetphi = jet->get_phi();
1138     _recojetmass = jet->get_mass();
1139     _recojetp = jet->get_p();
1140     _recojetenergy = jet->get_e();
1141     _recojetid = jet->get_id();
1142
1143     _truthjetid = 0;
1144     _truthjetp = 0;
1145     _truthjetphi = 0;
1146     _truthjeteta = 0;
1147     _truthjetpt = 0;
1148     _truthjetenergy = 0;
1149     _truthjetpx = 0;
1150     _truthjetpy = 0;
1151     _truthjetpz = 0;
1152
1153     //try to match reco jet with closest truth jet
1154     //this is A+A so we know the jet eval stack doesn't work
1155     //so just match the truth-reco jet pair by distance
1156     float closestjet = 9999;
1157     for (JetMap::Iter iter = truthjets->begin(); iter != truthjets->end(); ++iter)
1158     {
1159       Jet *truthjet = iter->second;
1160
1161       float thisjetpt = truthjet->get_pt();
1162       if (thisjetpt < minjetpt)
1163         continue;
1164
1165       float thisjeteta = truthjet->get_eta();
1166       float thisjetphi = truthjet->get_phi();
1167
1168       float dphi = recojetphi - thisjetphi;
1169       if (dphi > 3. * pi / 2.)
1170         dphi -= pi * 2.;
1171       if (dphi < -1. * pi / 2.)
1172         dphi += pi * 2.;
1173
1174       float deta = recojeteta - thisjeteta;
1175
1176       float dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
1177
1178       if (dR < recojets->get_par() && dR < closestjet)
1179       {
1180         _truthjetid = -9999;  //indicates matched with distance, not evalstack
1181         _truthjetp = truthjet->get_p();
1182         _truthjetphi = truthjet->get_phi();
1183         _truthjeteta = truthjet->get_eta();
1184         _truthjetpt = truthjet->get_pt();
1185         _truthjetenergy = truthjet->get_e();
1186         _truthjetpx = truthjet->get_px();
1187         _truthjetpy = truthjet->get_py();
1188         _truthjetpz = truthjet->get_pz();
1189         _truthjetmass = truthjet->get_mass();
1190       }
1191
1192       if (dR < closestjet)
1193         closestjet = dR;
1194     }
1195
1196     jetdphi = trig_phi - _recojetphi;
1197     if (jetdphi < pi2)
1198       jetdphi += 2. * pi;
1199     if (jetdphi > threepi2)
1200       jetdphi -= 2. * pi;
1201
1202     if (fabs(jetdphi) < 0.05)
1203       //don't care about matching the reco photon with itself
1204       continue;
1205
1206     jetdeta = trig_eta - _recojeteta;
1207     jetpout = _recojetpt * TMath::Sin(jetdphi);
1208
1209     isophot_jet_tree->Fill();
1210   }
1211 }
1212 void PhotonJet::GetRecoHadronsAndJets(RawCluster *trig,
1213                                       SvtxTrackMap *tracks,
1214                                       JetMap *jets,
1215                                       JetMap *trackedjets,
1216                                       JetRecoEval *recoeval,
1217                                       SvtxTrackEval *trackeval,
1218                                       PHG4TruthInfoContainer *alltruth,
1219                                       JetTruthEval *trutheval,
1220                                       JetMap *truthjets,
1221                                       GlobalVertex *vtx)
1222 {
1223   CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
1224   CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*trig, vertex);
1225
1226   float trig_phi = E_vec_cluster.getPhi();
1227   float trig_eta = E_vec_cluster.pseudoRapidity();
1228
1229   PHG4TruthInfoContainer::Range range = alltruth->GetPrimaryParticleRange();
1230
1231   //Match the reco photon with the truth photon
1232   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
1233        iter != range.second;
1234        ++iter)
1235   {
1236     PHG4Particle *truth = iter->second;
1237
1238     clustruthpid = truth->get_pid();
1239     if (clustruthpid == 22)
1240     {
1241       clustruthpx = truth->get_px();
1242       clustruthpy = truth->get_py();
1243       clustruthpz = truth->get_pz();
1244       clustruthenergy = truth->get_e();
1245       clustruthpt = sqrt(clustruthpx * clustruthpx + clustruthpy * clustruthpy);
1246       if (clustruthpt < mincluspt)
1247         continue;
1248
1249       TLorentzVector vec;
1250       vec.SetPxPyPzE(clustruthpx, clustruthpy, clustruthpz, clustruthenergy);
1251       clustruthphi = vec.Phi();
1252
1253       clustrutheta = TMath::ATanH(clustruthpz / clustruthenergy);
1254       //check for nans
1255       if (clustrutheta != clustrutheta)
1256         clustrutheta = -9999;
1257
1258       if (fabs(clustruthphi - trig_phi) > 0.02 || fabs(clustrutheta - trig_eta) > 0.02)
1259         continue;
1260
1261       //once the values are set, we've found the truth photon so just break out of this loop
1262       break;
1263     }
1264   }
1265
1266   if (eval_tracked_jets)
1267   {
1268     if (Verbosity() > 1)
1269     {
1270       cout << "evaluating tracked hadrons opposite the direct photon" << endl;
1271     }
1272     for (SvtxTrackMap::Iter iter = tracks->begin();
1273          iter != tracks->end();
1274          ++iter)
1275     {
1276       SvtxTrack *track = iter->second;
1277
1278       _tr_px = track->get_px();
1279       _tr_py = track->get_py();
1280       _tr_pz = track->get_pz();
1281       _tr_pt = sqrt(_tr_px * _tr_px + _tr_py * _tr_py);
1282       if (_tr_pt < 0.5)
1283         continue;
1284
1285       _tr_p = sqrt(_tr_px * _tr_px + _tr_py * _tr_py + _tr_pz * _tr_pz);
1286       _tr_phi = track->get_phi();
1287       _tr_eta = track->get_eta();
1288       _charge = track->get_charge();
1289       _chisq = track->get_chisq();
1290       //can find appropriate values to cut on for these later
1291       _ndf = track->get_ndf();
1292       _dca = track->get_dca();
1293
1294       _tr_x = track->get_x();
1295       _tr_y = track->get_y();
1296       _tr_z = track->get_z();
1297
1298       haddphi = trig_phi - _tr_phi;
1299       if (haddphi < pi2)
1300         haddphi += 2. * pi;
1301       if (haddphi > threepi2)
1302         haddphi -= 2. * pi;
1303
1304       //get the truth track info
1305       PHG4Particle *truthtrack = trackeval->max_truth_particle_by_nclusters(track);
1306       _truth_is_primary = alltruth->is_primary(truthtrack);
1307
1308       _truthtrackpx = truthtrack->get_px();
1309       _truthtrackpy = truthtrack->get_py();
1310       _truthtrackpz = truthtrack->get_pz();
1311       _truthtrackp = sqrt(truthtrackpx * truthtrackpx + truthtrackpy * truthtrackpy + truthtrackpz * truthtrackpz);
1312
1313       _truthtracke = truthtrack->get_e();
1314       TLorentzVector *Truthtrack = new TLorentzVector();
1315       Truthtrack->SetPxPyPzE(truthtrackpx, truthtrackpy, truthtrackpz, truthtracke);
1316       _truthtrackpt = Truthtrack->Pt();
1317       _truthtrackphi = Truthtrack->Phi();
1318       _truthtracketa = Truthtrack->Eta();
1319       _truthtrackpid = truthtrack->get_pid();
1320
1321       haddeta = trig_eta - _tr_eta;
1322
1323       hadpout = _tr_pt * TMath::Sin(haddphi);
1324
1325       isophot_had_tree->Fill();
1326     }
1327
1328     //now collect the away-sidet tracked jets from the direct photon
1329     for (JetMap::Iter iter = trackedjets->begin();
1330          iter != trackedjets->end();
1331          ++iter)
1332     {
1333       Jet *jet = iter->second;
1334       Jet *truthjet = recoeval->max_truth_jet_by_energy(jet);
1335
1336       _trecojetpt = jet->get_pt();
1337
1338       if (_trecojetpt < minjetpt)
1339         continue;
1340
1341       _trecojeteta = jet->get_eta();
1342       if (fabs(_trecojeteta) > 1.)
1343         continue;
1344
1345       _trecojetpx = jet->get_px();
1346       _trecojetpy = jet->get_py();
1347       _trecojetpz = jet->get_pz();
1348       _trecojetphi = jet->get_phi();
1349       _trecojetmass = jet->get_mass();
1350       _trecojetp = jet->get_p();
1351       _trecojetenergy = jet->get_e();
1352       _trecojetid = jet->get_id();
1353
1354       if (truthjet)
1355       {
1356         _ttruthjetid = truthjet->get_id();
1357         _ttruthjetp = truthjet->get_p();
1358         _ttruthjetphi = truthjet->get_phi();
1359         _ttruthjeteta = truthjet->get_eta();
1360         _ttruthjetpt = truthjet->get_pt();
1361         _ttruthjetenergy = truthjet->get_e();
1362         _ttruthjetpx = truthjet->get_px();
1363         _ttruthjetpy = truthjet->get_py();
1364         _ttruthjetpz = truthjet->get_pz();
1365       }
1366       //for some reason jeteval stack doesn't have matched jets, just set them to -999
1367       else
1368       {
1369         _ttruthjetid = -9999;
1370         _ttruthjetp = -9999;
1371         _ttruthjetphi = -9999;
1372         _ttruthjeteta = -9999;
1373         _ttruthjetpt = -9999;
1374         _ttruthjetenergy = -9999;
1375         _ttruthjetpx = -9999;
1376         _ttruthjetpy = -9999;
1377         _ttruthjetpz = -9999;
1378       }
1379
1380       tjetdphi = trig_phi - _trecojetphi;
1381       if (tjetdphi < pi2)
1382         tjetdphi += 2. * pi;
1383       if (tjetdphi > threepi2)
1384         tjetdphi -= 2. * pi;
1385
1386       if (fabs(tjetdphi) < 0.05)
1387         //just pairedthe photon with itself instead of with a jet, so skip it
1388         continue;
1389
1390       tjetdeta = trig_eta - _trecojeteta;
1391       tjetpout = _trecojetpt * TMath::Sin(jetdphi);
1392
1393       isophot_trackjet_tree->Fill();
1394     }
1395   }
1396
1397   //now collect awayside cluster jets
1398
1399   for (JetMap::Iter iter = jets->begin();
1400        iter != jets->end();
1401        ++iter)
1402   {
1403     Jet *jet = iter->second;
1404     Jet *truthjet = recoeval->max_truth_jet_by_energy(jet);
1405
1406     _recojetpt = jet->get_pt();
1407     if (_recojetpt < minjetpt)
1408       continue;
1409
1410     _recojeteta = jet->get_eta();
1411     if (_recojeteta < etalow || _recojeteta > etahigh)
1412       continue;
1413
1414     _recojetpx = jet->get_px();
1415     _recojetpy = jet->get_py();
1416     _recojetpz = jet->get_pz();
1417     _recojetphi = jet->get_phi();
1418     _recojetmass = jet->get_mass();
1419     _recojetp = jet->get_p();
1420     _recojetenergy = jet->get_e();
1421     _recojetid = jet->get_id();
1422     int pair_ntruthconstituents = 0;
1423     if (truthjet)
1424     {
1425       _truthjetid = truthjet->get_id();
1426       _truthjetp = truthjet->get_p();
1427       _truthjetphi = truthjet->get_phi();
1428       _truthjeteta = truthjet->get_eta();
1429       _truthjetpt = truthjet->get_pt();
1430       _truthjetenergy = truthjet->get_e();
1431       _truthjetpx = truthjet->get_px();
1432       _truthjetpy = truthjet->get_py();
1433       _truthjetpz = truthjet->get_pz();
1434       _truthjetmass = truthjet->get_mass();
1435
1436       //check that the jet isn't just a photon or something like this
1437       //needs at least 2 constituents and that 90% of jet isn't from one photon
1438       std::set<PHG4Particle *> truthjetcomp =
1439           trutheval->all_truth_particles(truthjet);
1440
1441       float truthjetenergysum = 0;
1442       float truthjethighestphoton = 0;
1443       for (std::set<PHG4Particle *>::iterator iter2 = truthjetcomp.begin();
1444            iter2 != truthjetcomp.end();
1445            ++iter2)
1446       {
1447         PHG4Particle *truthpart = *iter2;
1448         if (!truthpart)
1449         {
1450           cout << "no truth particles in the jet??" << endl;
1451           break;
1452         }
1453         pair_ntruthconstituents++;
1454
1455         int constpid = truthpart->get_pid();
1456         float conste = truthpart->get_e();
1457
1458         truthjetenergysum += conste;
1459
1460         if (constpid == 22)
1461         {
1462           if (conste > truthjethighestphoton)
1463             truthjethighestphoton = conste;
1464         }
1465       }
1466     }
1467
1468     else
1469     {
1470       _truthjetid = 0;
1471       _truthjetp = 0;
1472       _truthjetphi = 0;
1473       _truthjeteta = 0;
1474       _truthjetpt = 0;
1475       _truthjetenergy = 0;
1476       _truthjetpx = 0;
1477       _truthjetpy = 0;
1478       _truthjetpz = 0;
1479
1480       //if the jetevalstack can't find a good truth jet
1481       //try to match reco jet with closest truth jet
1482       float closestjet = 9999;
1483       for (JetMap::Iter iter = truthjets->begin(); iter != truthjets->end(); ++iter)
1484       {
1485         Jet *trutherjet = iter->second;
1486
1487         float thisjetpt = trutherjet->get_pt();
1488         if (thisjetpt < minjetpt)
1489           continue;
1490
1491         float thisjeteta = trutherjet->get_eta();
1492         float thisjetphi = trutherjet->get_phi();
1493
1494         float dphi = recojetphi - thisjetphi;
1495         if (dphi > 3. * pi / 2.)
1496           dphi -= pi * 2.;
1497         if (dphi < -1. * pi / 2.)
1498           dphi += pi * 2.;
1499
1500         float deta = recojeteta - thisjeteta;
1501
1502         float dR = sqrt(pow(dphi, 2.) + pow(deta, 2.));
1503
1504         if (dR < jets->get_par() && dR < closestjet)
1505         {
1506           _truthjetid = -9999;  //indicates matched with distance, not evalstack
1507           _truthjetp = trutherjet->get_p();
1508           _truthjetphi = trutherjet->get_phi();
1509           _truthjeteta = trutherjet->get_eta();
1510           _truthjetpt = trutherjet->get_pt();
1511           _truthjetenergy = trutherjet->get_e();
1512           _truthjetpx = trutherjet->get_px();
1513           _truthjetpy = trutherjet->get_py();
1514           _truthjetpz = trutherjet->get_pz();
1515           _truthjetmass = trutherjet->get_mass();
1516         }
1517
1518         if (dR < closestjet)
1519           closestjet = dR;
1520       }
1521     }
1522
1523     //make sure the jet has at least 3 constiutents
1524     if (!is_AA && pair_ntruthconstituents < 3)
1525       continue;
1526
1527     jetdphi = trig_phi - _recojetphi;
1528     if (jetdphi < pi2)
1529       jetdphi += 2. * pi;
1530     if (jetdphi > threepi2)
1531       jetdphi -= 2. * pi;
1532
1533     //just pairedthe photon with itself instead of with a jet, so skip it
1534     if (fabs(jetdphi) < 0.05)
1535       continue;
1536
1537     jetdeta = trig_eta - _recojeteta;
1538     jetpout = _recojetpt * TMath::Sin(jetdphi);
1539
1540     isophot_jet_tree->Fill();
1541   }
1542 }
1543
1544 float PhotonJet::ConeSum(RawCluster *cluster,
1545                          RawClusterContainer *cluster_container,
1546                          SvtxTrackMap *trackmap,
1547                          float coneradius,
1548                          GlobalVertex *vtx)
1549 {
1550   float energyptsum = 0;
1551
1552   CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
1553   CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
1554
1555   RawClusterContainer::ConstRange begin_end = cluster_container->getClusters();
1556   RawClusterContainer::ConstIterator clusiter;
1557
1558   for (clusiter = begin_end.first;
1559        clusiter != begin_end.second;
1560        ++clusiter)
1561   {
1562     RawCluster *conecluster = clusiter->second;
1563
1564     //check to make sure that the candidate isolated photon isn't being counted in the energy sum
1565     if (conecluster->get_energy() == cluster->get_energy())
1566       if (conecluster->get_phi() == cluster->get_phi())
1567         if (conecluster->get_z() == cluster->get_z())
1568           continue;
1569
1570     CLHEP::Hep3Vector E_vec_conecluster = RawClusterUtility::GetECoreVec(*conecluster, vertex);
1571
1572     float cone_pt = E_vec_conecluster.perp();
1573
1574     if (cone_pt < 0.2)
1575       continue;
1576
1577     float cone_e = conecluster->get_energy();
1578     float cone_eta = E_vec_conecluster.pseudoRapidity();
1579     float cone_phi = E_vec_conecluster.getPhi();
1580
1581     float dphi = cluster->get_phi() - cone_phi;
1582     if (dphi < -1 * pi)
1583       dphi += 2. * pi;
1584     if (dphi > pi)
1585       dphi -= 2. * pi;
1586
1587     float deta = E_vec_cluster.pseudoRapidity() - cone_eta;
1588
1589     float radius = sqrt(dphi * dphi + deta * deta);
1590
1591     if (radius < coneradius)
1592     {
1593       energyptsum += cone_e;
1594     }
1595   }
1596
1597   for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
1598   {
1599     SvtxTrack *track = iter->second;
1600
1601     float trackpx = track->get_px();
1602     float trackpy = track->get_py();
1603     float trackpt = sqrt(trackpx * trackpx + trackpy * trackpy);
1604     if (trackpt < 0.2)
1605       continue;
1606     float trackphi = track->get_phi();
1607     float tracketa = track->get_eta();
1608     float dphi = E_vec_cluster.getPhi() - trackphi;
1609     float deta = E_vec_cluster.pseudoRapidity() - tracketa;
1610     float radius = sqrt(dphi * dphi + deta * deta);
1611
1612     if (radius < coneradius)
1613     {
1614       energyptsum += trackpt;
1615     }
1616   }
1617
1618   return energyptsum;
1619 }
1620
1621 void PhotonJet::Set_Tree_Branches()
1622 {
1623   cluster_tree = new TTree("clustertree", "A tree with EMCal cluster information");
1624   cluster_tree->Branch("clus_energy", &clus_energy, "clus_energy/F");
1625   cluster_tree->Branch("clus_eta", &clus_eta, "clus_eta/F");
1626   cluster_tree->Branch("clus_phi", &clus_phi, "clus_phi/F");
1627   cluster_tree->Branch("clus_pt", &clus_pt, "clus_pt/F");
1628   cluster_tree->Branch("clus_theta", &clus_theta, "clus_theta/F");
1629   cluster_tree->Branch("clus_x", &clus_x, "clus_x/F");
1630   cluster_tree->Branch("clus_y", &clus_y, "clus_y/F");
1631   cluster_tree->Branch("clus_z", &clus_z, "clus_z/F");
1632   cluster_tree->Branch("clus_t", &clus_t, "clus_t/F");
1633   cluster_tree->Branch("clus_px", &clus_px, "clus_px/F");
1634   cluster_tree->Branch("clus_py", &clus_py, "clus_py/F");
1635   cluster_tree->Branch("clus_pz", &clus_pz, "clus_pz/F");
1636   cluster_tree->Branch("nevents", &nevents, "nevents/I");
1637   cluster_tree->Branch("process_id", &process_id, "process_id/I");
1638   cluster_tree->Branch("x1", &x1, "x1/F");
1639   cluster_tree->Branch("x2", &x2, "x2/F");
1640   cluster_tree->Branch("partid1", &partid1, "partid1/I");
1641   cluster_tree->Branch("partid2", &partid2, "partid2/I");
1642   cluster_tree->Branch("E_4x4", &E_4x4, "E_4x4/F");
1643   cluster_tree->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
1644   cluster_tree->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
1645   cluster_tree->Branch("E_2x2", &E_2x2, "E_2x2/F");
1646   cluster_tree->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
1647   cluster_tree->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
1648
1649   isolated_clusters = new TTree("isolated_clusters", "a tree with isolated EMCal clusters");
1650   isolated_clusters->Branch("clus_energy", &clus_energy, "clus_energy/F");
1651   isolated_clusters->Branch("clus_eta", &clus_eta, "clus_eta/F");
1652   isolated_clusters->Branch("clus_phi", &clus_phi, "clus_phi/F");
1653   isolated_clusters->Branch("clus_pt", &clus_pt, "clus_pt/F");
1654   isolated_clusters->Branch("clus_theta", &clus_theta, "clus_theta/F");
1655   isolated_clusters->Branch("clus_x", &clus_x, "clus_x/F");
1656   isolated_clusters->Branch("clus_y", &clus_y, "clus_y/F");
1657   isolated_clusters->Branch("clus_z", &clus_z, "clus_z/F");
1658   isolated_clusters->Branch("clus_t", &clus_t, "clus_t/F");
1659   isolated_clusters->Branch("clus_px", &clus_px, "clus_px/F");
1660   isolated_clusters->Branch("clus_py", &clus_py, "clus_py/F");
1661   isolated_clusters->Branch("clus_pz", &clus_pz, "clus_pz/F");
1662   isolated_clusters->Branch("nevents", &nevents, "nenvents/I");
1663   isolated_clusters->Branch("clustruthenergy", &clustruthenergy, "clustruthenergy/F");
1664   isolated_clusters->Branch("clustruthpt", &clustruthpt, "clustruthpt/F");
1665   isolated_clusters->Branch("clustruthphi", &clustruthphi, "clustruthphi/F");
1666   isolated_clusters->Branch("clustrutheta", &clustrutheta, "clustrutheta/F");
1667   isolated_clusters->Branch("clustruthpx", &clustruthpx, "clustruthpx/F");
1668   isolated_clusters->Branch("clustruthpy", &clustruthpy, "clustruthpy/F");
1669   isolated_clusters->Branch("clustruthpz", &clustruthpz, "clustruthpz/F");
1670   isolated_clusters->Branch("process_id", &process_id, "process_id/I");
1671   isolated_clusters->Branch("x1", &x1, "x1/F");
1672   isolated_clusters->Branch("x2", &x2, "x2/F");
1673   isolated_clusters->Branch("partid1", &partid1, "partid1/I");
1674   isolated_clusters->Branch("partid2", &partid2, "partid2/I");
1675   isolated_clusters->Branch("E_4x4", &E_4x4, "E_4x4/F");
1676   isolated_clusters->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
1677   isolated_clusters->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
1678   isolated_clusters->Branch("E_2x2", &E_2x2, "E_2x2/F");
1679   isolated_clusters->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
1680   isolated_clusters->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
1681
1682   tracktree = new TTree("tracktree", "a tree with svtx tracks");
1683   tracktree->Branch("tr_px", &tr_px, "tr_px/F");
1684   tracktree->Branch("tr_py", &tr_py, "tr_py/F");
1685   tracktree->Branch("tr_pz", &tr_pz, "tr_pz/F");
1686   tracktree->Branch("tr_pt", &tr_pt, "tr_pt/F");
1687   tracktree->Branch("tr_phi", &tr_phi, "tr_phi/F");
1688   tracktree->Branch("tr_eta", &tr_eta, "tr_eta/F");
1689   tracktree->Branch("charge", &charge, "charge/I");
1690   tracktree->Branch("chisq", &chisq, "chisq/F");
1691   tracktree->Branch("ndf", &ndf, "ndf/I");
1692   tracktree->Branch("dca", &dca, "dca/F");
1693   tracktree->Branch("tr_x", &tr_x, "tr_x/F");
1694   tracktree->Branch("tr_y", &tr_y, "tr_y/F");
1695   tracktree->Branch("tr_z", &tr_z, "tr_z/F");
1696   tracktree->Branch("nevents", &nevents, "nevents/i");
1697
1698   tracktree->Branch("truthtrackpx", &truthtrackpx, "truthtrackpx/F");
1699   tracktree->Branch("truthtrackpy", &truthtrackpy, "truthtrackpy/F");
1700   tracktree->Branch("truthtrackpz", &truthtrackpz, "truthtrackpz/F");
1701   tracktree->Branch("truthtrackp", &truthtrackp, "truthtrackp/F");
1702   tracktree->Branch("truthtracke", &truthtracke, "truthtracke/F");
1703   tracktree->Branch("truthtrackpt", &truthtrackpt, "truthtrackpt/F");
1704   tracktree->Branch("truthtrackphi", &truthtrackphi, "truthtrackphi/F");
1705   tracktree->Branch("truthtracketa", &truthtracketa, "truthtracketa/F");
1706   tracktree->Branch("truthtrackpid", &truthtrackpid, "truthtrackpid/I");
1707   tracktree->Branch("truth_is_primary", &truth_is_primary, "truth_is_primary/B");
1708   tracktree->Branch("process_id", &process_id, "process_id/I");
1709   tracktree->Branch("x1", &x1, "x1/F");
1710   tracktree->Branch("x2", &x2, "x2/F");
1711   tracktree->Branch("partid1", &partid1, "partid1/I");
1712   tracktree->Branch("partid2", &partid2, "partid2/I");
1713
1714   truthjettree = new TTree("truthjettree", "a tree with truth jets");
1715   truthjettree->Branch("ntruthconstituents", &ntruthconstituents, "ntruthconstituents/I");
1716   truthjettree->Branch("truthjetpt", &truthjetpt, "truthjetpt/F");
1717   truthjettree->Branch("truthjetpx", &truthjetpx, "truthjetpx/F");
1718   truthjettree->Branch("truthjetpy", &truthjetpy, "truthjetpy/F");
1719   truthjettree->Branch("truthjetpz", &truthjetpz, "truthjetpz/F");
1720   truthjettree->Branch("truthjetphi", &truthjetphi, "truthjetphi/F");
1721   truthjettree->Branch("truthjeteta", &truthjeteta, "truthjeteta/F");
1722   truthjettree->Branch("truthjetmass", &truthjetmass, "truthjetmass/F");
1723   truthjettree->Branch("truthjetp", &truthjetp, "truthjetp/F");
1724   truthjettree->Branch("truthjetenergy", &truthjetenergy, "truthjetenergy/F");
1725   truthjettree->Branch("nevents", &nevents, "nevents/I");
1726   truthjettree->Branch("process_id", &process_id, "process_id/I");
1727   truthjettree->Branch("x1", &x1, "x1/F");
1728   truthjettree->Branch("x2", &x2, "x2/F");
1729   truthjettree->Branch("partid1", &partid1, "partid1/I");
1730   truthjettree->Branch("partid2", &partid2, "partid2/I");
1731   truthjettree->Branch("E_4x4", &E_4x4, "E_4x4/F");
1732   truthjettree->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
1733   truthjettree->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
1734   truthjettree->Branch("E_2x2", &E_2x2, "E_2x2/F");
1735   truthjettree->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
1736   truthjettree->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
1737
1738   recojettree = new TTree("recojettree", "a tree with reco jets");
1739   recojettree->Branch("dR", &dR, "dR/F");
1740   recojettree->Branch("recojetpt", &recojetpt, "recojetpt/F");
1741   recojettree->Branch("recojetpx", &recojetpx, "recojetpx/F");
1742   recojettree->Branch("recojetpy", &recojetpy, "recojetpy/F");
1743   recojettree->Branch("recojetpz", &recojetpz, "recojetpz/F");
1744   recojettree->Branch("recojetphi", &recojetphi, "recojetphi/F");
1745   recojettree->Branch("recojeteta", &recojeteta, "recojeteta/F");
1746   recojettree->Branch("recojetmass", &recojetmass, "recojetmass/F");
1747   recojettree->Branch("recojetp", &recojetp, "recojetp/F");
1748   recojettree->Branch("recojetenergy", &recojetenergy, "recojetenergy/F");
1749   recojettree->Branch("nevents", &nevents, "nevents/I");
1750   recojettree->Branch("recojetid", &recojetid, "recojetid/F");
1751   recojettree->Branch("truthjetid", &truthjetid, "truthjetid/F");
1752   recojettree->Branch("truthjetp", &truthjetp, "truthjetp/F");
1753   recojettree->Branch("truthjetphi", &truthjetphi, "truthjetphi/F");
1754   recojettree->Branch("truthjeteta", &truthjeteta, "truthjeteta/F");
1755   recojettree->Branch("truthjetpt", &truthjetpt, "truthjetpt/F");
1756   recojettree->Branch("truthjetenergy", &truthjetenergy, "truthjetenergy/F");
1757   recojettree->Branch("truthjetpx", &truthjetpx, "truthjetpx/F");
1758   recojettree->Branch("truthjetpy", &truthjetpy, "truthjetpy/F");
1759   recojettree->Branch("truthjetpz", &truthjetpz, "truthjetpz/F");
1760   recojettree->Branch("process_id", &process_id, "process_id/I");
1761   recojettree->Branch("truthjetmass", &truthjetmass, "truthjetmass/F");
1762   recojettree->Branch("x1", &x1, "x1/F");
1763   recojettree->Branch("x2", &x2, "x2/F");
1764   recojettree->Branch("partid1", &partid1, "partid1/I");
1765   recojettree->Branch("partid2", &partid2, "partid2/I");
1766   recojettree->Branch("E_4x4", &E_4x4, "E_4x4/F");
1767   recojettree->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
1768   recojettree->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
1769   recojettree->Branch("E_2x2", &E_2x2, "E_2x2/F");
1770   recojettree->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
1771   recojettree->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
1772   recojettree->Branch("constituent_dphis","std::vector<float>",&constituent_dphis);
1773   recojettree->Branch("constituent_detas","std::vector<float>",&constituent_detas);
1774
1775
1776   isophot_jet_tree = new TTree("isophoton-jets",
1777                                "a tree with correlated isolated photons and jets");
1778   isophot_jet_tree->Branch("clus_energy", &clus_energy, "clus_energy/F");
1779   isophot_jet_tree->Branch("clus_eta", &clus_eta, "clus_eta/F");
1780   isophot_jet_tree->Branch("clus_phi", &clus_phi, "clus_phi/F");
1781   isophot_jet_tree->Branch("clus_pt", &clus_pt, "clus_pt/F");
1782   isophot_jet_tree->Branch("clus_theta", &clus_theta, "clus_theta/F");
1783   isophot_jet_tree->Branch("clus_x", &clus_x, "clus_x/F");
1784   isophot_jet_tree->Branch("clus_y", &clus_y, "clus_y/F");
1785   isophot_jet_tree->Branch("clus_z", &clus_z, "clus_z/F");
1786   isophot_jet_tree->Branch("clus_t", &clus_t, "clus_t/F");
1787   isophot_jet_tree->Branch("clus_px", &clus_px, "clus_px/F");
1788   isophot_jet_tree->Branch("clus_py", &clus_py, "clus_py/F");
1789   isophot_jet_tree->Branch("clus_pz", &clus_pz, "clus_pz/F");
1790   isophot_jet_tree->Branch("clustruthenergy", &clustruthenergy, "clustruthenergy/F");
1791   isophot_jet_tree->Branch("clustruthpt", &clustruthpt, "clustruthpt/F");
1792   isophot_jet_tree->Branch("clustruthphi", &clustruthphi, "clustruthphi/F");
1793   isophot_jet_tree->Branch("clustrutheta", &clustrutheta, "clustrutheta/F");
1794   isophot_jet_tree->Branch("clustruthpx", &clustruthpx, "clustruthpx/F");
1795   isophot_jet_tree->Branch("clustruthpy", &clustruthpy, "clustruthpy/F");
1796   isophot_jet_tree->Branch("clustruthpz", &clustruthpz, "clustruthpz/F");
1797   isophot_jet_tree->Branch("_recojetpt", &_recojetpt, "_recojetpt/F");
1798   isophot_jet_tree->Branch("_recojetpx", &_recojetpx, "_recojetpx/F");
1799   isophot_jet_tree->Branch("_recojetpy", &_recojetpy, "_recojetpy/F");
1800   isophot_jet_tree->Branch("_recojetpz", &_recojetpz, "_recojetpz/F");
1801   isophot_jet_tree->Branch("_recojetphi", &_recojetphi, "_recojetphi/F");
1802   isophot_jet_tree->Branch("_recojeteta", &_recojeteta, "_recojeteta/F");
1803   isophot_jet_tree->Branch("_recojetmass", &_recojetmass, "_recojetmass/F");
1804   isophot_jet_tree->Branch("_recojetp", &_recojetp, "_recojetp/F");
1805   isophot_jet_tree->Branch("_recojetenergy", &_recojetenergy, "_recojetenergy/F");
1806   isophot_jet_tree->Branch("jetdphi", &jetdphi, "jetdphi/F");
1807   isophot_jet_tree->Branch("jetdeta", &jetdeta, "jetdeta/F");
1808   isophot_jet_tree->Branch("jetpout", &jetpout, "jetpout/F");
1809   isophot_jet_tree->Branch("nevents", &nevents, "nevents/I");
1810   isophot_jet_tree->Branch("_recojetid", &_recojetid, "_recojetid/F");
1811   isophot_jet_tree->Branch("_truthjetid", &_truthjetid, "_truthjetid/F");
1812   isophot_jet_tree->Branch("_truthjetp", &_truthjetp, "_truthjetp/F");
1813   isophot_jet_tree->Branch("_truthjetphi", &_truthjetphi, "_truthjetphi/F");
1814   isophot_jet_tree->Branch("_truthjeteta", &_truthjeteta, "_truthjeteta/F");
1815   isophot_jet_tree->Branch("_truthjetpt", &_truthjetpt, "_truthjetpt/F");
1816   isophot_jet_tree->Branch("_truthjetenergy", &_truthjetenergy, "_truthjetenergy/F");
1817   isophot_jet_tree->Branch("_truthjetpx", &_truthjetpx, "_truthjetpx/F");
1818   isophot_jet_tree->Branch("_truthjetpy", &_truthjetpy, "_truthjetpy/F");
1819   isophot_jet_tree->Branch("_truthjetpz", &_truthjetpz, "_truthjetpz/F");
1820   isophot_jet_tree->Branch("process_id", &process_id, "process_id/I");
1821   isophot_jet_tree->Branch("x1", &x1, "x1/F");
1822   isophot_jet_tree->Branch("x2", &x2, "x2/F");
1823   isophot_jet_tree->Branch("partid1", &partid1, "partid1/I");
1824   isophot_jet_tree->Branch("partid2", &partid2, "partid2/I");
1825   isophot_jet_tree->Branch("E_4x4", &E_4x4, "E_4x4/F");
1826   isophot_jet_tree->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
1827   isophot_jet_tree->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
1828   isophot_jet_tree->Branch("E_2x2", &E_2x2, "E_2x2/F");
1829   isophot_jet_tree->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
1830   isophot_jet_tree->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
1831
1832   isophot_trackjet_tree = new TTree("isophoton-trackjets",
1833                                     "a tree with correlated isolated photons and track jets");
1834   isophot_trackjet_tree->Branch("clus_energy", &clus_energy, "clus_energy/F");
1835   isophot_trackjet_tree->Branch("clus_eta", &clus_eta, "clus_eta/F");
1836   isophot_trackjet_tree->Branch("clus_phi", &clus_phi, "clus_phi/F");
1837   isophot_trackjet_tree->Branch("clus_pt", &clus_pt, "clus_pt/F");
1838   isophot_trackjet_tree->Branch("clus_theta", &clus_theta, "clus_theta/F");
1839   isophot_trackjet_tree->Branch("clus_x", &clus_x, "clus_x/F");
1840   isophot_trackjet_tree->Branch("clus_y", &clus_y, "clus_y/F");
1841   isophot_trackjet_tree->Branch("clus_z", &clus_z, "clus_z/F");
1842   isophot_trackjet_tree->Branch("clus_t", &clus_t, "clus_t/F");
1843   isophot_trackjet_tree->Branch("clus_px", &clus_px, "clus_px/F");
1844   isophot_trackjet_tree->Branch("clus_py", &clus_py, "clus_py/F");
1845   isophot_trackjet_tree->Branch("clus_pz", &clus_pz, "clus_pz/F");
1846   isophot_trackjet_tree->Branch("clustruthenergy", &clustruthenergy, "clustruthenergy/F");
1847   isophot_trackjet_tree->Branch("clustruthpt", &clustruthpt, "clustruthpt/F");
1848   isophot_trackjet_tree->Branch("clustruthphi", &clustruthphi, "clustruthphi/F");
1849   isophot_trackjet_tree->Branch("clustrutheta", &clustrutheta, "clustrutheta/F");
1850   isophot_trackjet_tree->Branch("clustruthpx", &clustruthpx, "clustruthpx/F");
1851   isophot_trackjet_tree->Branch("clustruthpy", &clustruthpy, "clustruthpy/F");
1852   isophot_trackjet_tree->Branch("clustruthpz", &clustruthpz, "clustruthpz/F");
1853   isophot_trackjet_tree->Branch("_trecojetpt", &_trecojetpt, "_trecojetpt/F");
1854   isophot_trackjet_tree->Branch("_trecojetpx", &_trecojetpx, "_trecojetpx/F");
1855   isophot_trackjet_tree->Branch("_trecojetpy", &_trecojetpy, "_trecojetpy/F");
1856   isophot_trackjet_tree->Branch("_trecojetpz", &_trecojetpz, "_trecojetpz/F");
1857   isophot_trackjet_tree->Branch("_trecojetphi", &_trecojetphi, "_trecojetphi/F");
1858   isophot_trackjet_tree->Branch("_trecojeteta", &_trecojeteta, "_trecojeteta/F");
1859   isophot_trackjet_tree->Branch("_trecojetmass", &_trecojetmass, "_trecojetmass/F");
1860   isophot_trackjet_tree->Branch("_trecojetp", &_trecojetp, "_trecojetp/F");
1861   isophot_trackjet_tree->Branch("_trecojetenergy", &_trecojetenergy, "_trecojetenergy/F");
1862   isophot_trackjet_tree->Branch("tjetdphi", &tjetdphi, "tjetdphi/F");
1863   isophot_trackjet_tree->Branch("tjetdeta", &tjetdeta, "tjetdeta/F");
1864   isophot_trackjet_tree->Branch("tjetpout", &tjetpout, "tjetpout/F");
1865   isophot_trackjet_tree->Branch("nevents", &nevents, "nevents/I");
1866   isophot_trackjet_tree->Branch("_trecojetid", &_trecojetid, "_trecojetid/F");
1867   isophot_trackjet_tree->Branch("_ttruthjetid", &_ttruthjetid, "_ttruthjetid/F");
1868   isophot_trackjet_tree->Branch("_ttruthjetp", &_ttruthjetp, "_ttruthjetp/F");
1869   isophot_trackjet_tree->Branch("_ttruthjetphi", &_ttruthjetphi, "_ttruthjetphi/F");
1870   isophot_trackjet_tree->Branch("_ttruthjeteta", &_ttruthjeteta, "_ttruthjeteta/F");
1871   isophot_trackjet_tree->Branch("_ttruthjetpt", &_ttruthjetpt, "_ttruthjetpt/F");
1872   isophot_trackjet_tree->Branch("_ttruthjetenergy", &_ttruthjetenergy, "_ttruthjetenergy/F");
1873   isophot_trackjet_tree->Branch("_ttruthjetpx", &_ttruthjetpx, "_ttruthjetpx/F");
1874   isophot_trackjet_tree->Branch("_ttruthjetpy", &_ttruthjetpy, "_ttruthjetpy/F");
1875   isophot_trackjet_tree->Branch("_ttruthjetpz", &_ttruthjetpz, "_ttruthjetpz/F");
1876   isophot_trackjet_tree->Branch("process_id", &process_id, "process_id/I");
1877   isophot_trackjet_tree->Branch("x1", &x1, "x1/F");
1878   isophot_trackjet_tree->Branch("x2", &x2, "x2/F");
1879   isophot_trackjet_tree->Branch("partid1", &partid1, "partid1/I");
1880   isophot_trackjet_tree->Branch("partid2", &partid2, "partid2/I");
1881   isophot_trackjet_tree->Branch("E_4x4", &E_4x4, "E_4x4/F");
1882   isophot_trackjet_tree->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
1883   isophot_trackjet_tree->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
1884   isophot_trackjet_tree->Branch("E_2x2", &E_2x2, "E_2x2/F");
1885   isophot_trackjet_tree->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
1886   isophot_trackjet_tree->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
1887
1888   isophot_had_tree = new TTree("isophoton-hads", "a tree with correlated isolated photons and hadrons");
1889   isophot_had_tree->Branch("clus_energy", &clus_energy, "clus_energy/F");
1890   isophot_had_tree->Branch("clus_eta", &clus_eta, "clus_eta/F");
1891   isophot_had_tree->Branch("clus_phi", &clus_phi, "clus_phi/F");
1892   isophot_had_tree->Branch("clus_pt", &clus_pt, "clus_pt/F");
1893   isophot_had_tree->Branch("clus_theta", &clus_theta, "clus_theta/F");
1894   isophot_had_tree->Branch("clus_x", &clus_x, "clus_x/F");
1895   isophot_had_tree->Branch("clus_y", &clus_y, "clus_y/F");
1896   isophot_had_tree->Branch("clus_z", &clus_z, "clus_z/F");
1897   isophot_had_tree->Branch("clus_t", &clus_t, "clus_t/F");
1898   isophot_had_tree->Branch("clus_px", &clus_px, "clus_px/F");
1899   isophot_had_tree->Branch("clus_py", &clus_py, "clus_py/F");
1900   isophot_had_tree->Branch("clus_pz", &clus_pz, "clus_pz/F");
1901   isophot_had_tree->Branch("clustruthenergy", &clustruthenergy, "clustruthenergy/F");
1902   isophot_had_tree->Branch("clustruthpt", &clustruthpt, "clustruthpt/F");
1903   isophot_had_tree->Branch("clustruthphi", &clustruthphi, "clustruthphi/F");
1904   isophot_had_tree->Branch("clustrutheta", &clustrutheta, "clustrutheta/F");
1905   isophot_had_tree->Branch("clustruthpx", &clustruthpx, "clustruthpx/F");
1906   isophot_had_tree->Branch("clustruthpy", &clustruthpy, "clustruthpy/F");
1907   isophot_had_tree->Branch("clustruthpz", &clustruthpz, "clustruthpz/F");
1908
1909   isophot_had_tree->Branch("_tr_px", &_tr_px, "_tr_px/F");
1910   isophot_had_tree->Branch("_tr_py", &_tr_py, "_tr_py/F");
1911   isophot_had_tree->Branch("_tr_pz", &_tr_pz, "_tr_pz/F");
1912   isophot_had_tree->Branch("_tr_pt", &_tr_pt, "_tr_pt/F");
1913   isophot_had_tree->Branch("_tr_phi", &_tr_phi, "_tr_phi/F");
1914   isophot_had_tree->Branch("_tr_eta", &_tr_eta, "_tr_eta/F");
1915   isophot_had_tree->Branch("_charge", &_charge, "_charge/I");
1916   isophot_had_tree->Branch("_chisq", &_chisq, "_chisq/F");
1917   isophot_had_tree->Branch("_ndf", &_ndf, "_ndf/I");
1918   isophot_had_tree->Branch("_dca", &_dca, "_dca/F");
1919   isophot_had_tree->Branch("_tr_x", &_tr_x, "_tr_x/F");
1920   isophot_had_tree->Branch("_tr_y", &_tr_y, "_tr_y/F");
1921   isophot_had_tree->Branch("_tr_z", &_tr_z, "_tr_z/F");
1922   isophot_had_tree->Branch("haddphi", &haddphi, "haddphi/F");
1923   isophot_had_tree->Branch("hadpout", &hadpout, "hadpout/F");
1924   isophot_had_tree->Branch("haddeta", &haddeta, "haddeta/F");
1925   isophot_had_tree->Branch("nevents", &nevents, "nevents/I");
1926   isophot_had_tree->Branch("_truthtrackpx", &_truthtrackpx, "_truthtrackpx/F");
1927   isophot_had_tree->Branch("_truthtrackpy", &_truthtrackpy, "_truthtrackpy/F");
1928   isophot_had_tree->Branch("_truthtrackpz", &_truthtrackpz, "_truthtrackpz/F");
1929   isophot_had_tree->Branch("_truthtrackp", &_truthtrackp, "_truthtrackp/F");
1930   isophot_had_tree->Branch("_truthtracke", &_truthtracke, "_truthtracke/F");
1931   isophot_had_tree->Branch("_truthtrackpt", &_truthtrackpt, "_truthtrackpt/F");
1932   isophot_had_tree->Branch("_truthtrackphi", &_truthtrackphi, "_truthtrackphi/F");
1933   isophot_had_tree->Branch("_truthtracketa", &_truthtracketa, "_truthtracketa/F");
1934   isophot_had_tree->Branch("_truthtrackpid", &_truthtrackpid, "_truthtrackpid/I");
1935   isophot_had_tree->Branch("_truth_is_primary", &_truth_is_primary, "_truth_is_primary/B");
1936   isophot_had_tree->Branch("process_id", &process_id, "process_id/I");
1937   isophot_had_tree->Branch("x1", &x1, "x1/F");
1938   isophot_had_tree->Branch("x2", &x2, "x2/F");
1939   isophot_had_tree->Branch("partid1", &partid1, "partid1/I");
1940   isophot_had_tree->Branch("partid2", &partid2, "partid2/I");
1941   isophot_had_tree->Branch("E_4x4", &E_4x4, "E_4x4/F");
1942   isophot_had_tree->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
1943   isophot_had_tree->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
1944   isophot_had_tree->Branch("E_2x2", &E_2x2, "E_2x2/F");
1945   isophot_had_tree->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
1946   isophot_had_tree->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
1947
1948   truth_g4particles = new TTree("truthtree_g4", "a tree with all truth g4 particles");
1949   truth_g4particles->Branch("truthpx", &truthpx, "truthpx/F");
1950   truth_g4particles->Branch("truthpy", &truthpy, "truthpy/F");
1951   truth_g4particles->Branch("truthpz", &truthpz, "truthpz/F");
1952   truth_g4particles->Branch("truthp", &truthp, "truthp/F");
1953   truth_g4particles->Branch("truthenergy", &truthenergy, "truthenergy/F");
1954   truth_g4particles->Branch("truthphi", &truthphi, "truthphi/F");
1955   truth_g4particles->Branch("trutheta", &trutheta, "trutheta/F");
1956   truth_g4particles->Branch("truthpt", &truthpt, "truthpt/F");
1957   truth_g4particles->Branch("truthpid", &truthpid, "truthpid/I");
1958   truth_g4particles->Branch("nevents", &nevents, "nevents/I");
1959   truth_g4particles->Branch("process_id", &process_id, "process_id/I");
1960   truth_g4particles->Branch("x1", &x1, "x1/F");
1961   truth_g4particles->Branch("x2", &x2, "x2/F");
1962   truth_g4particles->Branch("partid1", &partid1, "partid1/I");
1963   truth_g4particles->Branch("partid2", &partid2, "partid2/I");
1964   truth_g4particles->Branch("E_4x4", &E_4x4, "E_4x4/F");
1965   truth_g4particles->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
1966   truth_g4particles->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
1967   truth_g4particles->Branch("E_2x2", &E_2x2, "E_2x2/F");
1968   truth_g4particles->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
1969   truth_g4particles->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
1970
1971   truthtree = new TTree("truthtree", "a tree with all truth pythia particles");
1972   truthtree->Branch("truthpx", &truthpx, "truthpx/F");
1973   truthtree->Branch("truthpy", &truthpy, "truthpy/F");
1974   truthtree->Branch("truthpz", &truthpz, "truthpz/F");
1975   truthtree->Branch("truthp", &truthp, "truthp/F");
1976   truthtree->Branch("truthenergy", &truthenergy, "truthenergy/F");
1977   truthtree->Branch("truthphi", &truthphi, "truthphi/F");
1978   truthtree->Branch("trutheta", &trutheta, "trutheta/F");
1979   truthtree->Branch("truthpt", &truthpt, "truthpt/F");
1980   truthtree->Branch("truthpid", &truthpid, "truthpid/I");
1981   truthtree->Branch("nevents", &nevents, "nevents/I");
1982   truthtree->Branch("numparticlesinevent", &numparticlesinevent, "numparticlesinevent/I");
1983   truthtree->Branch("process_id", &process_id, "process_id/I");
1984   truthtree->Branch("x1", &x1, "x1/F");
1985   truthtree->Branch("x2", &x2, "x2/F");
1986   truthtree->Branch("partid1", &partid1, "partid1/I");
1987   truthtree->Branch("partid2", &partid2, "partid2/I");
1988
1989   event_tree = new TTree("eventtree", "A tree with 2 to 2 event info");
1990   event_tree->Branch("mpi", &mpi, "mpi/F");
1991   event_tree->Branch("x1", &x1, "x1/F");
1992   event_tree->Branch("x2", &x2, "x2/F");
1993   event_tree->Branch("partid1", &partid1, "partid1/I");
1994   event_tree->Branch("partid2", &partid2, "partid2/I");
1995   event_tree->Branch("process_id", &process_id, "process_id/I");
1996   event_tree->Branch("nevents", &nevents, "nevents/I");
1997   event_tree->Branch("cluseventenergy", &cluseventenergy, "cluseventenergy/F");
1998   event_tree->Branch("cluseventeta", &cluseventeta, "cluseventeta/F");
1999   event_tree->Branch("cluseventpt", &cluseventpt, "cluseventpt/F");
2000   event_tree->Branch("cluseventphi", &cluseventphi, "cluseventphi/F");
2001   event_tree->Branch("hardest_jetpt", &hardest_jetpt, "hardest_jetpt/F");
2002   event_tree->Branch("hardest_jetphi", &hardest_jetphi, "hardest_jetphi/F");
2003   event_tree->Branch("hardest_jeteta", &hardest_jeteta, "hardest_jeteta/F");
2004   event_tree->Branch("hardest_jetenergy", &hardest_jetenergy, "hardest_jetenergy/F");
2005   event_tree->Branch("E_4x4", &E_4x4, "E_4x4/F");
2006   event_tree->Branch("phi_4x4", &phi_4x4, "phi_4x4/F");
2007   event_tree->Branch("eta_4x4", &eta_4x4, "eta_4x4/F");
2008   event_tree->Branch("E_2x2", &E_2x2, "E_2x2/F");
2009   event_tree->Branch("phi_2x2", &phi_2x2, "phi_2x2/F");
2010   event_tree->Branch("eta_2x2", &eta_2x2, "eta_2x2/F");
2011 }
2012
2013 void PhotonJet::initialize_values()
2014 {
2015   E_4x4 = 0;
2016   phi_4x4 = 0;
2017   eta_4x4 = 0;
2018   E_2x2 = 0;
2019   phi_2x2 = 0;
2020   eta_2x2 = 0;
2021
2022   event_tree = 0;
2023   tree = 0;
2024   cluster_tree = 0;
2025   truth_g4particles = 0;
2026   truthtree = 0;
2027   isolated_clusters = 0;
2028   tracktree = 0;
2029   truthjettree = 0;
2030   recojettree = 0;
2031   isophot_jet_tree = 0;
2032   isophot_trackjet_tree = 0;
2033   isophot_had_tree = 0;
2034
2035   histo = 0;
2036   dR = 0;
2037   percent_photon_h = 0;
2038   _truthjetmass = -999;
2039   clus_energy = -999;
2040   clus_eta = -999;
2041   clus_phi = -999;
2042   clus_pt = -999;
2043   clus_px = -999;
2044   clus_py = -999;
2045   clus_pz = -999;
2046   clus_theta = -999;
2047   clus_x = -999;
2048   clus_y = -999;
2049   clus_z = -999;
2050   clus_t = -999;
2051
2052   tr_px = -999;
2053   tr_py = -999;
2054   tr_pz = -999;
2055   tr_p = -999;
2056   tr_pt = -999;
2057   tr_phi = -999;
2058   tr_eta = -999;
2059   charge = -999;
2060   chisq = -999;
2061   ndf = -999;
2062   dca = -999;
2063   tr_x = -999;
2064   tr_y = -999;
2065   tr_z = -999;
2066   truthtrackpx = -999;
2067   truthtrackpy = -999;
2068   truthtrackpz = -999;
2069   truthtrackp = -999;
2070   truthtracke = -999;
2071   truthtrackpt = -999;
2072   truthtrackphi = -999;
2073   truthtracketa = -999;
2074   truthtrackpid = -999;
2075   truth_is_primary = false;
2076
2077   truthjetpt = -999;
2078   truthjetpx = -999;
2079   truthjetpy = -999;
2080   truthjetpz = -999;
2081   truthjetphi = -999;
2082   truthjeteta = -999;
2083   truthjetmass = -999;
2084   truthjetp = -999;
2085   truthjetenergy = -999;
2086   recojetpt = -999;
2087   recojetpx = -999;
2088   recojetpy = -999;
2089   recojetpz = -999;
2090   recojetphi = -999;
2091   recojeteta = -999;
2092   recojetmass = -999;
2093   recojetp = -999;
2094   recojetenergy = -999;
2095   recojetid = -999;
2096   truthjetid = -999;
2097
2098   _recojetid = -999;
2099   _recojetpt = -999;
2100   _recojetpx = -999;
2101   _recojetpy = -999;
2102   _recojetpz = -999;
2103   _recojetphi = -999;
2104   _recojeteta = -999;
2105   _recojetmass = -999;
2106   _recojetp = -999;
2107   _recojetenergy = -999;
2108   jetdphi = -999;
2109   jetpout = -999;
2110   jetdeta = -999;
2111   _truthjetid = -999;
2112   _truthjetpt = -999;
2113   _truthjetpx = -999;
2114   _truthjetpy = -999;
2115   _truthjetpz = -999;
2116   _truthjetphi = -999;
2117   _truthjeteta = -999;
2118   _truthjetmass = -999;
2119   _truthjetp = -999;
2120   _truthjetenergy = -999;
2121
2122   _trecojetid = -999;
2123   _trecojetpt = -999;
2124   _trecojetpx = -999;
2125   _trecojetpy = -999;
2126   _trecojetpz = -999;
2127   _trecojetphi = -999;
2128   _trecojeteta = -999;
2129   _trecojetmass = -999;
2130   _trecojetp = -999;
2131   _trecojetenergy = -999;
2132   tjetdphi = -999;
2133   tjetpout = -999;
2134   tjetdeta = -999;
2135   _ttruthjetid = -999;
2136   _ttruthjetpt = -999;
2137   _ttruthjetpx = -999;
2138   _ttruthjetpy = -999;
2139   _ttruthjetpz = -999;
2140   _ttruthjetphi = -999;
2141   _ttruthjeteta = -999;
2142   _ttruthjetmass = -999;
2143   _ttruthjetp = -999;
2144   _ttruthjetenergy = -999;
2145
2146   _tr_px = -999;
2147   _tr_py = -999;
2148   _tr_pz = -999;
2149   _tr_p = -999;
2150   _tr_pt = -999;
2151   _tr_phi = -999;
2152   _tr_eta = -999;
2153   _charge = -999;
2154   _chisq = -999;
2155   _ndf = -999;
2156   _dca = -999;
2157   _tr_x = -999;
2158   _tr_y = -999;
2159   _tr_z = -999;
2160   haddphi = -999;
2161   hadpout = -999;
2162   haddeta = -999;
2163   _truth_is_primary = false;
2164   _truthtrackpx = -999;
2165   _truthtrackpy = -999;
2166   _truthtrackpz = -999;
2167   _truthtrackp = -999;
2168   _truthtracke = -999;
2169   _truthtrackpt = -999;
2170   _truthtrackphi = -999;
2171   _truthtracketa = -999;
2172   _truthtrackpid = -999;
2173
2174   truthpx = -999;
2175   truthpy = -999;
2176   truthpz = -999;
2177   truthp = -999;
2178   truthphi = -999;
2179   trutheta = -999;
2180   truthpt = -999;
2181   truthenergy = -999;
2182   truthpid = -999;
2183   numparticlesinevent = -999;
2184   process_id = -999;
2185
2186   clustruthpx = -999;
2187   clustruthpy = -999;
2188   clustruthpz = -999;
2189   clustruthenergy = -999;
2190   clustruthpt = -999;
2191   clustruthphi = -999;
2192   clustrutheta = -999;
2193   clustruthpid = -999;
2194
2195   beam_energy = 0;
2196   x1 = 0;
2197   x2 = 0;
2198   partid1 = 0;
2199   partid2 = 0;
2200
2201   hardest_jetpt = 0;
2202   hardest_jetphi = 0;
2203   hardest_jeteta = 0;
2204   hardest_jetenergy = 0;
2205 }