PhotonJet/Photons/Photons.C

0001 //This is an analysis package designed to collect single photons
0002 //They can be embedded or not, and can be in either the central
0003 //or forward arms - there are trees for each
0004
0005
0006 #include "Photons.h"
0007
0008 #include <fun4all/Fun4AllServer.h>
0009 #include <g4main/PHG4Hit.h>
0010 #include <g4main/PHG4Particle.h>
0011 #include <g4main/PHG4TruthInfoContainer.h>
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/getClass.h>
0014
0015 #include <TLorentzVector.h>
0016 #include <g4jets/Jet.h>
0017 #include <g4jets/JetMap.h>
0018 #include <iostream>
0019
0020 #include <calobase/RawCluster.h>
0021 #include <calobase/RawClusterContainer.h>
0022 #include <calobase/RawClusterUtility.h>
0023 #include <calobase/RawTower.h>
0024 #include <calobase/RawTowerContainer.h>
0025 #include <calobase/RawTowerGeom.h>
0026 #include <calobase/RawTowerGeomContainer.h>
0027 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
0028 #include <g4eval/JetEvalStack.h>
0029
0030 #include <g4vertex/GlobalVertex.h>
0031 #include <g4vertex/GlobalVertexMap.h>
0032
0033 #include <g4eval/SvtxEvalStack.h>
0034 #include <sstream>
0035
0036 #include <HepMC/GenEvent.h>
0037 #include <HepMC/GenVertex.h>
0038 #include <phhepmc/PHHepMCGenEvent.h>
0039 using namespace std;
0040
0041 #include <cassert>
0042 #include <iostream>
0043
0044 Photons::Photons(const std::string &name)
0045   : SubsysReco("PHOTONS")
0046 {
0047   outfilename = name;
0048   //initialize global variables to -999 so that they have a spot in memory
0049   initialize_to_zero();
0050
0051   //add other initializers here
0052
0053
0054   //default central arm
0055   _etalow = -1;
0056   _etahi = 1;
0057
0058   nevents = 0;
0059
0060   //default no hijing embedding
0061   _embed = 0;
0062 }
0063
0064 int Photons::Init(PHCompositeNode *topnode)
0065 {
0066   file = new TFile(outfilename.c_str(), "RECREATE");
0067
0068   histo = new TH1F("histo", "histo", 100, -3, 3);
0069
0070   tree = new TTree("tree", "a tree");
0071   tree->Branch("nevents", &nevents, "nevents/I");
0072
0073   //set all the tree branches
0074   Set_Tree_Branches();
0075
0076   return 0;
0077 }
0078
0079 int Photons::process_event(PHCompositeNode *topnode)
0080 {
0081   if (nevents % 10 == 0)
0082     cout << "at event number " << nevents << endl;
0083
0084
0085
0086   //get the nodes from the NodeTree
0087   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topnode, "G4TruthInfo");
0088   //Raw clusters
0089   RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_CEMC");
0090
0091   //Recalib clusters
0092   RawClusterContainer *recal_clusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_POS_COR_CEMC");
0093
0094   RawClusterContainer *fclusters = 0;
0095   if (_etalow > 1)
0096   {
0097     fclusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_FEMC");
0098   }
0099
0100   RawClusterContainer *hcalin_clusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_HCALIN");
0101
0102   RawTowerContainer *_towers = findNode::getClass<RawTowerContainer>(topnode, "TOWER_CALIB_CEMC");
0103
0104   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topnode, "GlobalVertexMap");
0105   if (!vertexmap)
0106   {
0107     cout << "Photons::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;
0108     assert(vertexmap);  // force quit
0109
0110     return 0;
0111   }
0112
0113   if (vertexmap->empty())
0114   {
0115     cout << "Photons::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;
0116     return 0;
0117   }
0118
0119   GlobalVertex *vtx = vertexmap->begin()->second;
0120   if (vtx == nullptr) return 0;
0121
0122   if (!recal_clusters)
0123   {
0124     cout << "no recalibrated cemc clusters, bailing" << endl;
0125     return 0;
0126   }
0127
0128   if (!_towers)
0129   {
0130     cout << "NO CALIBRATED CEMC TOWERS, BAILING" << endl;
0131     return 0;
0132   }
0133   if (!fclusters && _etalow > 1)
0134   {
0135     cout << "not forward cluster info" << endl;
0136     return 0;
0137   }
0138   if (!truthinfo)
0139   {
0140     cout << "no truth track info" << endl;
0141     return 0;
0142   }
0143   if (!clusters)
0144   {
0145     cout << "no cluster info" << endl;
0146     return 0;
0147   }
0148
0149   if (!hcalin_clusters)
0150   {
0151     cout << "no hcal in cluster info" << endl;
0152     return 0;
0153   }
0154
0155   if (Verbosity() > 1)
0156   {
0157     cout << "Getting truth particles" << endl;
0158   }
0159
0160   PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0161
0162   for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0163   {
0164     PHG4Particle *truth = iter->second;
0165
0166     const int this_embed_id = truthinfo->isEmbeded(truth->get_track_id());
0167
0168     //check that the truth particle is from the event generator
0169     //and not from the hijing background.
0170     //if it is from the hijing background this_embed_id == 0
0171     if (this_embed_id != 1 && _embed)
0172       continue;
0173     truthpid = truth->get_pid();
0174     truthpx = truth->get_px();
0175     truthpy = truth->get_py();
0176     truthpz = truth->get_pz();
0177     truthp = sqrt(truthpx * truthpx + truthpy * truthpy + truthpz * truthpz);
0178     truthenergy = truth->get_e();
0179
0180     TLorentzVector vec;
0181     vec.SetPxPyPzE(truthpx, truthpy, truthpz, truthenergy);
0182
0183     truthpt = sqrt(truthpx * truthpx + truthpy * truthpy);
0184
0185     truthphi = vec.Phi();
0186     trutheta = vec.Eta();
0187
0188     truth_g4particles->Fill();
0189   }
0190
0191   /***********************************************
0192
0193   GET THE HCAL_INNER CLUSTERS
0194   (for looking for leakage between emcal towers and/or tunneling)
0195   This map not be necessary in the future if we don't have an
0196   inner HCal :/
0197   ************************************************/
0198
0199   RawClusterContainer::ConstRange begin_end_hcal = hcalin_clusters->getClusters();
0200   RawClusterContainer::ConstIterator hcaliter;
0201
0202   if (Verbosity() > 1)
0203   {
0204     cout << "Getting inner HCal clusters for energy leakage studies" << endl;
0205   }
0206
0207   for (hcaliter = begin_end_hcal.first;
0208        hcaliter != begin_end_hcal.second;
0209        ++hcaliter)
0210   {
0211     RawCluster *cluster = hcaliter->second;
0212     CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0213     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0214     hcal_energy = E_vec_cluster.mag();
0215     hcal_eta = E_vec_cluster.pseudoRapidity();
0216     hcal_theta = E_vec_cluster.getTheta();
0217     hcal_pt = E_vec_cluster.perp();
0218     hcal_phi = E_vec_cluster.getPhi();
0219
0220     if (hcal_pt < 0.5)
0221       continue;
0222
0223     TLorentzVector *clus = new TLorentzVector();
0224     clus->SetPtEtaPhiE(hcal_pt, hcal_eta, hcal_phi, hcal_energy);
0225     float dumarray[4];
0226     clus->GetXYZT(dumarray);
0227     hcal_x = dumarray[0];
0228     hcal_y = dumarray[1];
0229     hcal_z = dumarray[2];
0230     hcal_t = dumarray[3];
0231
0232     hcal_px = hcal_pt * TMath::Cos(hcal_phi);
0233     hcal_py = hcal_pt * TMath::Sin(hcal_phi);
0234     hcal_pz = sqrt(hcal_energy * hcal_energy - hcal_px * hcal_px * hcal_py * hcal_py);
0235
0236     //find the associated truth high pT photon with this reconstructed
0237     //hcal cluster
0238
0239     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0240     {
0241       PHG4Particle *truth = iter->second;
0242
0243       hclustruthpid = truth->get_pid();
0244       if (hclustruthpid == 22)
0245       {
0246         hclustruthpx = truth->get_px();
0247         hclustruthpy = truth->get_py();
0248         hclustruthpz = truth->get_pz();
0249         hclustruthenergy = truth->get_e();
0250         hclustruthpt = sqrt(clustruthpx * clustruthpx + clustruthpy * clustruthpy);
0251         if (hclustruthpt < 0.5)
0252           continue;
0253
0254         TLorentzVector vec;
0255         vec.SetPxPyPzE(hclustruthpx, hclustruthpy, hclustruthpz, hclustruthenergy);
0256         hclustruthphi = vec.Phi();
0257         hclustrutheta = vec.Eta();
0258         //once found it break out
0259         break;
0260       }
0261     }
0262
0263     inhcal_tree->Fill();
0264   }
0265
0266   /***********************************************
0267
0268   GET THE FORWARD EMCAL CLUSTERS
0269
0270   ************************************************/
0271
0272   if (_etalow > 1)
0273   {
0274     if (Verbosity() > 1)
0275     {
0276       cout << "getting forward EMCal clusters" << endl;
0277     }
0278
0279     RawClusterContainer::ConstRange fclus = fclusters->getClusters();
0280     RawClusterContainer::ConstIterator fclusiter;
0281
0282     for (fclusiter = fclus.first;
0283      fclusiter != fclus.second;
0284      ++fclusiter)
0285     {
0286       RawCluster *cluster = fclusiter->second;
0287
0288       CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0289       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0290       fclusenergy = E_vec_cluster.mag();
0291       fclus_eta = E_vec_cluster.pseudoRapidity();
0292       fclus_theta = E_vec_cluster.getTheta();
0293       fclus_pt = E_vec_cluster.perp();
0294       fclus_phi = E_vec_cluster.getPhi();
0295
0296       fclus_px = fclus_pt * TMath::Cos(fclus_phi);
0297       fclus_py = fclus_pt * TMath::Sin(fclus_phi);
0298       fclus_pz = fclusenergy * TMath::Cos(fclus_theta);
0299
0300
0301
0302       //find the truth photon that corresponds to this event
0303       for (PHG4TruthInfoContainer::ConstIterator fiter = range.first;
0304        fiter != range.second;
0305        ++fiter)
0306       {
0307         PHG4Particle *truth = fiter->second;
0308
0309         fclustruthpid = truth->get_pid();
0310     //can run for photons or electrons
0311         if (fclustruthpid == 22 || abs(fclustruthpid) == 11)
0312         {
0313           fclustruthpx = truth->get_px();
0314           fclustruthpy = truth->get_py();
0315           fclustruthpz = truth->get_pz();
0316           fclustruthenergy = truth->get_e();
0317           fclustruthpt = sqrt(fclustruthpx * fclustruthpx + fclustruthpy * fclustruthpy);
0318           if (fclustruthpt < 0.3)
0319             continue;
0320
0321           TLorentzVector vec;
0322           vec.SetPxPyPzE(fclustruthpx, fclustruthpy, fclustruthpz, fclustruthenergy);
0323           fclustruthphi = vec.Phi();
0324           fclustrutheta = vec.Eta();
0325           //once found it break out
0326           break;
0327         }
0328       }
0329
0330       fcluster_tree->Fill();
0331     }
0332   }
0333
0334
0335
0336
0337
0338
0339   /***********************************************
0340
0341   GET THE POSITION RECALIBRATED EMCAL CLUSTERS
0342
0343   ************************************************/
0344
0345   RawClusterContainer::ConstRange rbegin_end = recal_clusters->getClusters();
0346   RawClusterContainer::ConstIterator rclusiter;
0347
0348   if (Verbosity() > 1)
0349     cout << "Getting the position recalibrated clusters" << endl;
0350
0351   //check that we are analyzing central arms
0352   //position correction doesn't exist for forward arms
0353   if (_etahi < 1.1)
0354   {
0355     for (rclusiter = rbegin_end.first;
0356      rclusiter != rbegin_end.second;
0357      ++rclusiter)
0358       {
0359
0360       RawCluster *cluster = rclusiter->second;
0361
0362       CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0363       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0364       rclus_energy = E_vec_cluster.mag();
0365       rclus_eta = E_vec_cluster.pseudoRapidity();
0366       rclus_theta = E_vec_cluster.getTheta();
0367       rclus_pt = E_vec_cluster.perp();
0368       rclus_phi = E_vec_cluster.getPhi();
0369
0370       if (rclus_pt < mincluspt)
0371         continue;
0372       if(Verbosity() > 1)
0373     cout<<"passed recal pt cut"<<endl;
0374       TLorentzVector *clus = new TLorentzVector();
0375       clus->SetPtEtaPhiE(rclus_pt, rclus_eta, rclus_phi, rclus_energy);
0376
0377       float dumarray[4];
0378       clus->GetXYZT(dumarray);
0379       rclus_x = dumarray[0];
0380       rclus_y = dumarray[1];
0381       rclus_z = dumarray[2];
0382       rclus_t = dumarray[3];
0383
0384       rclus_px = rclus_pt * TMath::Cos(rclus_phi);
0385       rclus_py = rclus_pt * TMath::Sin(rclus_phi);
0386       rclus_pz = sqrt(rclus_energy * rclus_energy - rclus_px * rclus_px - rclus_py * rclus_py);
0387       rclus_ecore = cluster->get_ecore();
0388       rclus_prob = cluster->get_prob();
0389       rclus_chi2 = cluster->get_chi2();
0390
0391       //find the associated truth high pT photon with this reconstructed photon
0392
0393       for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0394       {
0395         PHG4Particle *truth = iter->second;
0396
0397         clustruthpid = truth->get_pid();
0398         const int this_embed_id = truthinfo->isEmbeded(truth->get_track_id());
0399
0400     //if it is embedded make sure that it is from the event
0401     //generation and not the HIJING background
0402     if (this_embed_id != 1 && _embed)
0403           continue;
0404
0405     //can run for photons or electrons in the EMCal
0406         if (clustruthpid == 22 || fabs(clustruthpid) == 11)
0407         {
0408           clustruthpx = truth->get_px();
0409           clustruthpy = truth->get_py();
0410           clustruthpz = truth->get_pz();
0411           clustruthenergy = truth->get_e();
0412           clustruthpt = sqrt(clustruthpx * clustruthpx + clustruthpy * clustruthpy);
0413
0414           if (clustruthpt < mincluspt)
0415             continue;
0416
0417           TLorentzVector vec;
0418           vec.SetPxPyPzE(clustruthpx, clustruthpy, clustruthpz, clustruthenergy);
0419           clustruthphi = vec.Phi();
0420           clustrutheta = vec.Eta();
0421           //once found it break out
0422       if(Verbosity() > 1)
0423         cout<<"found recal truth photon"<<endl;
0424           break;
0425         }
0426       }
0427       if(Verbosity() > 1)
0428     cout<<"filling recal cluster tree"<<endl;
0429
0430       recal_cluster_tree->Fill();
0431       }
0432   }
0433
0434
0435
0436
0437
0438   /***********************************************
0439
0440   GET THE REGULAR, NON POSITION CORRECTED EMCAL CLUSTERS
0441
0442   ************************************************/
0443
0444   RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0445   RawClusterContainer::ConstIterator clusiter;
0446
0447   if (Verbosity() > 1)
0448     cout << "Get the non-position recalibrated clusters" << endl;
0449
0450   for (clusiter = begin_end.first;
0451        clusiter != begin_end.second;
0452        ++clusiter)
0453   {
0454
0455     RawCluster *cluster = clusiter->second;
0456
0457     CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0458     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0459     clus_energy = E_vec_cluster.mag();
0460     clus_eta = E_vec_cluster.pseudoRapidity();
0461     clus_theta = E_vec_cluster.getTheta();
0462     clus_pt = E_vec_cluster.perp();
0463     clus_phi = E_vec_cluster.getPhi();
0464
0465     clus_ecore = cluster->get_ecore();
0466     clus_chi2 = cluster->get_chi2();
0467     clus_prob = cluster->get_prob();
0468
0469     if (clus_pt < mincluspt)
0470       continue;
0471
0472     if(Verbosity() > 1)
0473       cout<<"passed cluster pt cut"<<endl;
0474
0475     TLorentzVector *clus = new TLorentzVector();
0476     clus->SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clus_energy);
0477
0478     float dumarray[4];
0479     clus->GetXYZT(dumarray);
0480     clus_x = dumarray[0];
0481     clus_y = dumarray[1];
0482     clus_z = dumarray[2];
0483     clus_t = dumarray[3];
0484
0485     clus_px = clus_pt * TMath::Cos(clus_phi);
0486     clus_py = clus_pt * TMath::Sin(clus_phi);
0487     clus_pz = sqrt(clus_energy * clus_energy - clus_px * clus_px - clus_py * clus_py);
0488
0489     //find the associated truth high pT photon with this reconstructed photon
0490
0491     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0492     {
0493       PHG4Particle *truth = iter->second;
0494
0495       clustruthpid = truth->get_pid();
0496       const int this_embed_id = truthinfo->isEmbeded(truth->get_track_id());
0497       if (this_embed_id != 1 && _embed)
0498         continue;
0499       if (clustruthpid == 22 || fabs(clustruthpid) == 11)
0500       {
0501         clustruthpx = truth->get_px();
0502         clustruthpy = truth->get_py();
0503         clustruthpz = truth->get_pz();
0504         clustruthenergy = truth->get_e();
0505         clustruthpt = sqrt(clustruthpx * clustruthpx + clustruthpy * clustruthpy);
0506
0507         if (clustruthpt < mincluspt)
0508           continue;
0509
0510         TLorentzVector vec;
0511         vec.SetPxPyPzE(clustruthpx, clustruthpy, clustruthpz, clustruthenergy);
0512         clustruthphi = vec.Phi();
0513         clustrutheta = vec.Eta();
0514         //once found it break out
0515
0516         break;
0517       }
0518     }
0519
0520     /***********************************************
0521
0522    DO THE EMCAL RECALIBRATION (if need be, recals are in database as of January 2018)
0523    I'm leaving the code here commented out in case it ever needs to be used
0524    again for new recalibrations e.g. for a new clusterizer
0525
0526   ************************************************/
0527     /*
0528     std::vector<float> toweretas;
0529     std::vector<float> towerphis;
0530     std::vector<float> towerenergies;
0531
0532     RawCluster::TowerConstRange towerbegend = cluster->get_towers();
0533     RawCluster::TowerConstIterator toweriter;
0534
0535     //loop over the towers in the cluster
0536     for(toweriter=towerbegend.first;
0537     toweriter!=towerbegend.second;
0538     ++toweriter){
0539
0540       RawTower *tower = _towers->getTower(toweriter->first);
0541
0542       int towereta = tower->get_bineta();
0543       int towerphi = tower->get_binphi();
0544       double towerenergy = tower->get_energy();
0545
0546       //put the etabin, phibin, and energy into the corresponding vectors
0547       toweretas.push_back(towereta);
0548       towerphis.push_back(towerphi);
0549       towerenergies.push_back(towerenergy);
0550
0551     }
0552
0553     int ntowers = toweretas.size();
0554
0555     //loop over the towers to determine the energy
0556     //weighted eta and phi position of the cluster
0557
0558     float etamult=0;
0559     float etasum=0;
0560     float phimult=0;
0561     float phisum=0;
0562
0563     for(int j=0; j<ntowers; j++){
0564
0565       float energymult = towerenergies.at(j)*toweretas.at(j);
0566       etamult+=energymult;
0567       etasum+=towerenergies.at(j);
0568
0569       energymult = towerenergies.at(j)*towerphis.at(j);
0570       phimult+=energymult;
0571       phisum+=towerenergies.at(j);
0572     }
0573     float avgphi = phimult/phisum;
0574     float avgeta = etamult/etasum;
0575
0576     fmodphi = fmod(avgphi,2.);
0577     fmodeta = fmod(avgeta,2.);
0578
0579
0580
0581     cluster_tree->Fill();
0582
0583
0584   }
0585
0586
0587
0588
0589
0590   //TOWER GEOMETRY, if needed can use to check the actual limits of the towers
0591
0592
0593    RawTowerGeomContainer *_towergeoms = findNode::getClass<RawTowerGeomContainer>(topnode,"TOWERGEOM_CEMC");
0594
0595    if(!_towergeoms){
0596      cout<<"no tower geometry, bailing"<<endl;
0597      return 0;
0598    }
0599
0600    RawTowerGeomContainer::ConstRange towbegend =
0601      _towergeoms->get_tower_geometries();
0602    RawTowerGeomContainer::ConstIterator towiter;
0603    for(towiter=towbegend.first;
0604        towiter!=towbegend.second;
0605        ++towiter){
0606
0607
0608      RawTowerGeom *geom = towiter->second;
0609
0610
0611
0612      double eta = geom->get_eta();
0613      double phi = geom->get_phi();
0614      double center_x = geom->get_center_x();
0615      double center_y = geom->get_center_y();
0616      double center_z = geom->get_center_z();
0617      double size_x = geom->get_size_x();
0618      double size_y = geom->get_size_y();
0619      double size_z = geom->get_size_z();
0620
0621      cout<<eta<<"   "<<phi<<"   "<<center_x<<","<<center_y<<","<<center_z<<"     "<<size_x<<","<<size_y<<","<<size_z<<endl;
0622
0623     */
0624
0625
0626
0627     cluster_tree->Fill();
0628   }
0629
0630   if(Verbosity() > 1)
0631     cout<<"Finished event in Photons package"<<endl;
0632
0633   nevents++;
0634   tree->Fill();
0635   return 0;
0636 }
0637
0638 int Photons::End(PHCompositeNode *topnode)
0639 {
0640   std::cout << " DONE PROCESSING " << endl;
0641
0642   file->Write();
0643   file->Close();
0644   return 0;
0645 }
0646
0647 void Photons::Set_Tree_Branches()
0648 {
0649   inhcal_tree = new TTree("hcalclustree", "a tree with inner hcal cluster information");
0650   inhcal_tree->Branch("hcal_energy", &hcal_energy, "hcal_energy/F");
0651   inhcal_tree->Branch("hcal_eta", &hcal_eta, "hcal_eta/F");
0652   inhcal_tree->Branch("hcal_phi", &hcal_phi, "hcal_phi/F");
0653   inhcal_tree->Branch("hcal_pt", &hcal_pt, "hcal_pt/F");
0654   inhcal_tree->Branch("hcal_theta", &hcal_theta, "hcal_theta/F");
0655   inhcal_tree->Branch("hcal_x", &hcal_x, "hcal_x/F");
0656   inhcal_tree->Branch("hcal_y", &hcal_y, "hcal_y/F");
0657   inhcal_tree->Branch("hcal_z", &hcal_z, "hcal_z/F");
0658   inhcal_tree->Branch("hcal_t", &hcal_t, "hcal_t/F");
0659   inhcal_tree->Branch("hcal_px", &hcal_px, "hcal_px/F");
0660   inhcal_tree->Branch("hcal_py", &hcal_py, "hcal_py/F");
0661   inhcal_tree->Branch("hcal_pz", &hcal_pz, "hcal_pz/F");
0662   inhcal_tree->Branch("nevents", &nevents, "nevents/I");
0663   inhcal_tree->Branch("hclustruthpx", &hclustruthpx, "hclustruthpx/F");
0664   inhcal_tree->Branch("hclustruthpy", &hclustruthpy, "hclustruthpy/F");
0665   inhcal_tree->Branch("hclustruthpz", &hclustruthpz, "hclustruthpz/F");
0666   inhcal_tree->Branch("hclustruthenergy", &hclustruthenergy, "hclustruthenergy/F");
0667   inhcal_tree->Branch("hclustruthpt", &hclustruthpt, "hclustruthpt/F");
0668   inhcal_tree->Branch("hclustruthphi", &hclustruthphi, "hclustruthphi/F");
0669   inhcal_tree->Branch("hclustrutheta", &hclustrutheta, "hclustrutheta/F");
0670
0671   cluster_tree = new TTree("clustertree", "A tree with EMCal cluster information");
0672   cluster_tree->Branch("clus_energy", &clus_energy, "clus_energy/F");
0673   cluster_tree->Branch("clus_eta", &clus_eta, "clus_eta/F");
0674   cluster_tree->Branch("clus_phi", &clus_phi, "clus_phi/F");
0675   cluster_tree->Branch("clus_pt", &clus_pt, "clus_pt/F");
0676   cluster_tree->Branch("clus_theta", &clus_theta, "clus_theta/F");
0677   cluster_tree->Branch("clus_x", &clus_x, "clus_x/F");
0678   cluster_tree->Branch("clus_y", &clus_y, "clus_y/F");
0679   cluster_tree->Branch("clus_z", &clus_z, "clus_z/F");
0680   cluster_tree->Branch("clus_t", &clus_t, "clus_t/F");
0681   cluster_tree->Branch("clus_px", &clus_px, "clus_px/F");
0682   cluster_tree->Branch("clus_py", &clus_py, "clus_py/F");
0683   cluster_tree->Branch("clus_pz", &clus_pz, "clus_pz/F");
0684   cluster_tree->Branch("nevents", &nevents, "nevents/I");
0685   cluster_tree->Branch("clus_ecore", &clus_ecore, "clus_ecore/F");
0686   cluster_tree->Branch("clus_chi2", &clus_chi2, "clus_chi2/F");
0687   cluster_tree->Branch("clustruthpx", &clustruthpx, "clustruthpx/F");
0688   cluster_tree->Branch("clustruthpy", &clustruthpy, "clustruthpy/F");
0689   cluster_tree->Branch("clustruthpz", &clustruthpz, "clustruthpz/F");
0690   cluster_tree->Branch("clustruthenergy", &clustruthenergy, "clustruthenergy/F");
0691   cluster_tree->Branch("clustruthpt", &clustruthpt, "clustruthpt/F");
0692   cluster_tree->Branch("clustruthphi", &clustruthphi, "clustruthphi/F");
0693   cluster_tree->Branch("clustrutheta", &clustrutheta, "clustrutheta/F");
0694   cluster_tree->Branch("fmodphi", &fmodphi, "fmodphi/F");
0695   cluster_tree->Branch("fmodeta", &fmodeta, "fmodeta/F");
0696
0697   recal_cluster_tree = new TTree("recalclustertree", "A tree with EMCal recalib cluster information");
0698   recal_cluster_tree->Branch("rclus_energy", &rclus_energy, "rclus_energy/F");
0699   recal_cluster_tree->Branch("rclus_eta", &rclus_eta, "rclus_eta/F");
0700   recal_cluster_tree->Branch("rclus_phi", &rclus_phi, "rclus_phi/F");
0701   recal_cluster_tree->Branch("rclus_pt", &rclus_pt, "rclus_pt/F");
0702   recal_cluster_tree->Branch("rclus_theta", &rclus_theta, "rclus_theta/F");
0703   recal_cluster_tree->Branch("rclus_ecore", &rclus_ecore, "rclus_ecore/F");
0704   recal_cluster_tree->Branch("rclus_chi2", &rclus_chi2, "rclus_chi2/F");
0705   recal_cluster_tree->Branch("rclus_x", &rclus_x, "rclus_x/F");
0706   recal_cluster_tree->Branch("rclus_y", &rclus_y, "rclus_y/F");
0707   recal_cluster_tree->Branch("rclus_z", &rclus_z, "rclus_z/F");
0708   recal_cluster_tree->Branch("rclus_t", &rclus_t, "rclus_t/F");
0709   recal_cluster_tree->Branch("rclus_px", &rclus_px, "rclus_px/F");
0710   recal_cluster_tree->Branch("rclus_py", &rclus_py, "rclus_py/F");
0711   recal_cluster_tree->Branch("rclus_pz", &rclus_pz, "rclus_pz/F");
0712   recal_cluster_tree->Branch("nevents", &nevents, "nevents/I");
0713   recal_cluster_tree->Branch("clustruthpx", &clustruthpx, "clustruthpx/F");
0714   recal_cluster_tree->Branch("clustruthpy", &clustruthpy, "clustruthpy/F");
0715   recal_cluster_tree->Branch("clustruthpz", &clustruthpz, "clustruthpz/F");
0716   recal_cluster_tree->Branch("clustruthenergy", &clustruthenergy, "clustruthenergy/F");
0717   recal_cluster_tree->Branch("clustruthpt", &clustruthpt, "clustruthpt/F");
0718   recal_cluster_tree->Branch("clustruthphi", &clustruthphi, "clustruthphi/F");
0719   recal_cluster_tree->Branch("clustrutheta", &clustrutheta, "clustrutheta/F");
0720
0721   fcluster_tree = new TTree("fclustertree", "A tree with FEMCal cluster information");
0722   fcluster_tree->Branch("fclusenergy", &fclusenergy, "fclusenergy/F");
0723   fcluster_tree->Branch("fclus_eta", &fclus_eta, "fclus_eta/F");
0724   fcluster_tree->Branch("fclus_phi", &fclus_phi, "fclus_phi/F");
0725   fcluster_tree->Branch("fclus_pt", &fclus_pt, "fclus_pt/F");
0726   fcluster_tree->Branch("fclus_theta", &fclus_theta, "fclus_theta/F");
0727   fcluster_tree->Branch("fclus_px", &fclus_px, "fclus_px/F");
0728   fcluster_tree->Branch("fclus_py", &fclus_py, "fclus_py/F");
0729   fcluster_tree->Branch("fclus_pz", &fclus_pz, "fclus_pz/F");
0730   fcluster_tree->Branch("nevents", &nevents, "nevents/I");
0731   fcluster_tree->Branch("fclustruthpx", &fclustruthpx, "fclustruthpx/F");
0732   fcluster_tree->Branch("fclustruthpy", &fclustruthpy, "fclustruthpy/F");
0733   fcluster_tree->Branch("fclustruthpz", &fclustruthpz, "fclustruthpz/F");
0734   fcluster_tree->Branch("fclustruthenergy", &fclustruthenergy, "fclustruthenergy/F");
0735   fcluster_tree->Branch("fclustruthpt", &fclustruthpt, "fclustruthpt/F");
0736   fcluster_tree->Branch("fclustruthphi", &fclustruthphi, "fclustruthphi/F");
0737   fcluster_tree->Branch("fclustrutheta", &fclustrutheta, "fclustrutheta/F");
0738
0739   truth_g4particles = new TTree("truthtree_g4", "a tree with all truth g4 particles");
0740   truth_g4particles->Branch("truthpx", &truthpx, "truthpx/F");
0741   truth_g4particles->Branch("truthpy", &truthpy, "truthpy/F");
0742   truth_g4particles->Branch("truthpz", &truthpz, "truthpz/F");
0743   truth_g4particles->Branch("truthp", &truthp, "truthp/F");
0744   truth_g4particles->Branch("truthenergy", &truthenergy, "truthenergy/F");
0745   truth_g4particles->Branch("truthphi", &truthphi, "truthphi/F");
0746   truth_g4particles->Branch("trutheta", &trutheta, "trutheta/F");
0747   truth_g4particles->Branch("truthpt", &truthpt, "truthpt/F");
0748   truth_g4particles->Branch("truthpid", &truthpid, "truthpid/I");
0749   truth_g4particles->Branch("nevents", &nevents, "nevents/I");
0750 }
0751
0752 void Photons::initialize_to_zero()
0753 {
0754   file = 0;
0755   tree = 0;
0756   cluster_tree = 0;
0757   truth_g4particles = 0;
0758   inhcal_tree = 0;
0759   fcluster_tree = 0;
0760
0761   recal_cluster_tree = 0;
0762   rclus_energy = 0;
0763   rclus_eta = 0;
0764   rclus_phi = 0;
0765   rclus_pt = 0;
0766   rclus_px = 0;
0767   rclus_pz = 0;
0768   rclus_py = 0;
0769   rclus_theta = 0;
0770   rclus_x = 0;
0771   rclus_y = 0;
0772   rclus_z = 0;
0773   rclus_t = 0;
0774
0775   nevents = 0;
0776   histo = 0;
0777   hcal_energy = 0;
0778   hcal_eta = 0;
0779   hcal_phi = 0;
0780   hcal_pt = 0;
0781   hcal_px = 0;
0782   hcal_py = 0;
0783   hcal_pz = 0;
0784   hcal_theta = 0;
0785   hcal_x = 0;
0786   hcal_y = 0;
0787   hcal_z = 0;
0788   hcal_t = 0;
0789
0790   clus_energy = 0;
0791   clus_eta = 0;
0792   clus_phi = 0;
0793   clus_pt = 0;
0794   clus_px = 0;
0795   clus_py = 0;
0796   clus_pz = 0;
0797   clus_theta = 0;
0798   clus_x = 0;
0799   clus_y = 0;
0800   clus_z = 0;
0801   clus_t = 0;
0802   fmodphi = 0;
0803   fmodeta = 0;
0804
0805   truthpx = 0;
0806   truthpy = 0;
0807   truthpz = 0;
0808   truthp = 0;
0809   truthphi = 0;
0810   trutheta = 0;
0811   truthpt = 0;
0812   truthenergy = 0;
0813   truthpid = 0;
0814   numparticlesinevent = 0;
0815   process_id = 0;
0816
0817   clustruthpx = 0;
0818   clustruthpy = 0;
0819   clustruthpz = 0;
0820   clustruthenergy = 0;
0821   clustruthpt = 0;
0822   clustruthphi = 0;
0823   clustrutheta = 0;
0824   clustruthpid = 0;
0825   hclustruthpx = 0;
0826   hclustruthpy = 0;
0827   hclustruthpz = 0;
0828   hclustruthenergy = 0;
0829   hclustruthpt = 0;
0830   hclustruthphi = 0;
0831   hclustrutheta = 0;
0832   hclustruthpid = 0;
0833
0834   fclusenergy = 0;
0835   fclus_eta = 0;
0836   fclus_phi = 0;
0837   fclus_theta = 0;
0838   fclus_pt = 0;
0839   fclustruthpid = 0;
0840   fclustruthpx = 0;
0841   fclustruthpy = 0;
0842   fclustruthpz = 0;
0843   fclustruthenergy = 0;
0844   fclustruthpt = 0;
0845   fclustruthphi = 0;
0846   fclustrutheta = 0;
0847   fclus_px = 0;
0848   fclus_py = 0;
0849   fclus_pz = 0;
0850 }