PhotonJet/Photons/Forward_pi0s.C

0001 #include "Forward_pi0s.h"
0002
0003 #include <fun4all/Fun4AllServer.h>
0004 #include <g4main/PHG4Hit.h>
0005 #include <g4main/PHG4Particle.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/getClass.h>
0009
0010 #include <TLorentzVector.h>
0011 #include <calobase/RawTower.h>
0012 #include <calobase/RawTowerContainer.h>
0013 #include <calobase/RawTowerGeom.h>
0014 #include <calobase/RawTowerGeomContainer.h>
0015 #include <g4jets/Jet.h>
0016 #include <g4jets/JetMap.h>
0017 #include <iostream>
0018
0019 #include <calobase/RawCluster.h>
0020 #include <calobase/RawClusterContainer.h>
0021 #include <calobase/RawClusterUtility.h>
0022 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
0023 #include <g4eval/JetEvalStack.h>
0024 #include <g4eval/SvtxEvalStack.h>
0025
0026 #include <g4vertex/GlobalVertex.h>
0027 #include <g4vertex/GlobalVertexMap.h>
0028 #include <sstream>
0029
0030 #include <HepMC/GenEvent.h>
0031 #include <HepMC/GenVertex.h>
0032 #include <phhepmc/PHHepMCGenEvent.h>
0033 using namespace std;
0034
0035 #include <cassert>
0036 #include <iostream>
0037
0038 Forward_pi0s::Forward_pi0s(const std::string &name)
0039   : SubsysReco("FORWARD_PI0S")
0040 {
0041   outfilename = name;
0042
0043   //add other initializers here
0044   //default use isocone algorithm
0045   use_isocone = 1;
0046
0047   //default central arm
0048   _etalow = -1;
0049   _etahi = 1;
0050
0051   //default to only central arm
0052   _useforwardarm = 0;
0053
0054   mincluspt = 0.3;
0055
0056   nevents = 0;
0057   //default use 0.3 jet cone
0058   jet_cone_size = 0.3;
0059
0060 }
0061
0062 int Forward_pi0s::Init(PHCompositeNode *topnode)
0063 {
0064   file = new TFile(outfilename.c_str(), "RECREATE");
0065
0066   histo = new TH1F("histo", "histo", 100, -3, 3);
0067
0068   tree = new TTree("tree", "a tree");
0069   tree->Branch("nevents", &nevents, "nevents/I");
0070
0071   Set_Tree_Branches();
0072   initialize_things();
0073
0074   return 0;
0075 }
0076
0077 int Forward_pi0s::process_event(PHCompositeNode *topnode)
0078 {
0079   nevents++;
0080   cout << "at event number " << nevents << endl;
0081
0082   //reset the truth variables
0083   tphote2 = -999;
0084   tphotpx2 = -999;
0085   tphotpy2 = -999;
0086   tphotpz2 = -999;
0087   tphotpid2 = -999;
0088   tphotparentid2 = -999;
0089   tphotpt2 = -999;
0090   tphotphi2 = -999;
0091   tphoteta2 = -999;
0092
0093   tphote1 = -999;
0094   tphotpx1 = -999;
0095   tphotpy1 = -999;
0096   tphotpz1 = -999;
0097   tphotpid1 = -999;
0098   tphotparentid1 = -999;
0099   tphotpt1 = -999;
0100   tphotphi1 = -999;
0101   tphoteta1 = -999;
0102
0103   tpi0e = -999;
0104   tpi0px = -999;
0105   tpi0py = -999;
0106   tpi0pz = -999;
0107   tpi0pid = -999;
0108   tpi0pt = -999;
0109   tpi0phi = -999;
0110   tpi0eta = -999;
0111
0112   //reset cluster variables
0113   fclusenergy = -999;
0114   fclus_eta = -999;
0115   fclus_phi = -999;
0116   fclus_theta = -999;
0117   fclus_pt = -999;
0118   fclus_px = -999;
0119   fclus_py = -999;
0120   fclus_pz = -999;
0121
0122
0123   clus_energy = -999;
0124   clus_eta = -999;
0125   clus_theta = -999;
0126   clus_pt = -999;
0127   clus_phi = -999;
0128   clus_px = -999;
0129   clus_py = -999;
0130   clus_pz = -999;
0131
0132
0133
0134
0135
0136
0137   //get the nodes from the NodeTree
0138
0139   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topnode, "G4TruthInfo");
0140   RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_CEMC");
0141   RawClusterContainer *fclusters = 0;
0142   if (_useforwardarm)
0143   {
0144     fclusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_FEMC");
0145   }
0146   PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0147
0148   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topnode, "GlobalVertexMap");
0149   if (!vertexmap)
0150   {
0151     cout << "Forward_pi0s::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;
0152     assert(vertexmap);  // force quit
0153
0154     return 0;
0155   }
0156
0157   if (vertexmap->empty())
0158   {
0159     cout << "Forward_pi0s::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;
0160     return 0;
0161   }
0162
0163   GlobalVertex *vtx = vertexmap->begin()->second;
0164   if (vtx == nullptr) return 0;
0165
0166
0167   if (!fclusters && _useforwardarm)
0168   {
0169     cout << "not forward cluster info" << endl;
0170     return 0;
0171   }
0172   if (!truthinfo)
0173   {
0174     cout << "no truth track info" << endl;
0175     return 0;
0176   }
0177
0178   if (!clusters)
0179   {
0180     cout << "no cluster info" << endl;
0181     return 0;
0182   }
0183
0184
0185   for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0186   {
0187     PHG4Particle *truth = iter->second;
0188
0189     truthpx = truth->get_px();
0190     truthpy = truth->get_py();
0191     truthpz = truth->get_pz();
0192     truthp = sqrt(truthpx * truthpx + truthpy * truthpy + truthpz * truthpz);
0193     truthenergy = truth->get_e();
0194
0195     TLorentzVector vec;
0196     vec.SetPxPyPzE(truthpx, truthpy, truthpz, truthenergy);
0197
0198     truthpt = sqrt(truthpx * truthpx + truthpy * truthpy);
0199
0200     truthphi = vec.Phi();
0201     trutheta = vec.Eta();
0202     truthpid = truth->get_pid();
0203
0204     //get the pi0 information
0205     if (truthpid == 111)
0206     {
0207       tpi0e = truth->get_e();
0208       tpi0px = truth->get_px();
0209       tpi0py = truth->get_py();
0210       tpi0pz = truth->get_pz();
0211       tpi0pid = truth->get_pid();
0212       tpi0pt = sqrt(tpi0px * tpi0px + tpi0py * tpi0py);
0213       TLorentzVector pi0;
0214       pi0.SetPxPyPzE(tpi0px, tpi0py, tpi0pz, tpi0e);
0215       tpi0phi = pi0.Phi();
0216       tpi0eta = pi0.Eta();
0217     }
0218
0219     truth_g4particles->Fill();
0220   }
0221
0222   if (Verbosity() > 1)
0223     cout << "Truth pi0 energy is :" << tpi0e << endl;
0224
0225   int firstphotonflag = 0;
0226   int firstfirstphotonflag = 0;
0227   int secondphotonflag = 0;
0228   int secondsecondphotonflag = 0;
0229   PHG4TruthInfoContainer::Range ranger = truthinfo->GetSecondaryParticleRange();
0230   //This sets the global variables in the if statements to the pion first decay photon truth values and second decay photon truth values
0231   for (PHG4TruthInfoContainer::ConstIterator iter = ranger.first; iter != ranger.second; ++iter)
0232   {
0233     PHG4Particle *truth = iter->second;
0234
0235     int dumtruthpid = truth->get_pid();
0236     int dumparentid = truth->get_parent_id();
0237
0238
0239     //if the parent is the pi0 and its a photon and we haven't marked one yet
0240     if (dumparentid == 1 && dumtruthpid == 22 && !firstphotonflag)
0241     {
0242       tphote1 = truth->get_e();
0243       tphotpx1 = truth->get_px();
0244       tphotpy1 = truth->get_py();
0245       tphotpz1 = truth->get_pz();
0246       tphotpid1 = truth->get_pid();
0247       tphotparentid1 = truth->get_parent_id();
0248       tphotpt1 = sqrt(tphotpx1 * tphotpx1 + tphotpy1 * tphotpy1);
0249
0250       TLorentzVector vec;
0251       vec.SetPxPyPzE(tphotpx1, tphotpy1, tphotpz1, tphote1);
0252
0253       tphotphi1 = vec.Phi();
0254       tphoteta1 = vec.Eta();
0255       firstphotonflag = 1;
0256     }
0257
0258     //if the parent is the pi0 and its a photon and we have
0259     //marked the first photon but not the second photon
0260     if (dumparentid == 1 &&
0261         dumtruthpid == 22 &&
0262         firstphotonflag && firstfirstphotonflag)
0263     {
0264       tphote2 = truth->get_e();
0265       tphotpx2 = truth->get_px();
0266       tphotpy2 = truth->get_py();
0267       tphotpz2 = truth->get_pz();
0268       tphotpid2 = truth->get_pid();
0269       tphotparentid2 = truth->get_parent_id();
0270       tphotpt2 = sqrt(tphotpx2 * tphotpx2 + tphotpy2 * tphotpy2);
0271
0272       TLorentzVector vec;
0273       vec.SetPxPyPzE(tphotpx2, tphotpy2, tphotpz2, tphote2);
0274
0275       tphotphi2 = vec.Phi();
0276       tphoteta2 = vec.Eta();
0277       secondphotonflag = 1;
0278     }
0279
0280     //this is a decay to photon + ee pair i.e. dalitz decay, lets just skip it
0281     if (dumparentid == 1 &&
0282         firstphotonflag && secondphotonflag && secondsecondphotonflag)
0283     {
0284       cout << "I AM GOING TO SKIP THIS EVENT BECAUSE IT APPEARS TO BE A DALITZ DECAY" << endl;
0285       cout << "THERE WERE 3 PARTICLES MARKED WITH PARENT ID == 1" << endl;
0286
0287       return 0;
0288     }
0289
0290     //we need these extra flags because otherwise the dalitz check
0291     //will always have first and secondphoton flags marked as true
0292     //so this allows the dalitz check to occur after marking the
0293     //second truth photon as a decay component
0294     if (firstphotonflag)
0295       firstfirstphotonflag = 1;
0296     if (secondphotonflag)
0297       secondsecondphotonflag = 1;
0298   }
0299   if (Verbosity() > 1)
0300   {
0301     cout << "truth decay photon 1 energy | phi | eta " << tphote1
0302          << "  |   " << tphotphi1 << "   |   " << tphoteta1 << endl;
0303     cout << "truth decay photon 2 energy | phi | eta " << tphote2
0304          << "  |   " << tphotphi2 << "   |   " << tphoteta2 << endl;
0305   }
0306
0307
0308
0309
0310
0311
0312   /***********************************************
0313
0314   GET THE FORWARD EMCAL CLUSTERS
0315
0316   ************************************************/
0317
0318   if (_useforwardarm)
0319   {
0320     RawClusterContainer::ConstRange fclus = fclusters->getClusters();
0321     RawClusterContainer::ConstIterator fclusiter;
0322
0323     for (fclusiter = fclus.first; fclusiter != fclus.second; ++fclusiter)
0324     {
0325       RawCluster *cluster = fclusiter->second;
0326
0327       CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0328       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0329       fclusenergy = E_vec_cluster.mag();
0330       fclus_eta = E_vec_cluster.pseudoRapidity();
0331       fclus_theta = E_vec_cluster.getTheta();
0332       fclus_pt = E_vec_cluster.perp();
0333       fclus_phi = E_vec_cluster.getPhi();
0334
0335
0336       if (fclusenergy < mincluspt)
0337         continue;
0338
0339
0340       fclus_px = fclus_pt * TMath::Cos(fclus_phi);
0341       fclus_py = fclus_pt * TMath::Sin(fclus_phi);
0342       fclus_pz = fclusenergy * TMath::Cos(fclus_theta);
0343
0344       //reset second cluster variables
0345       fclusenergy2 = -999;
0346       fclus_eta2 = -999;
0347       fclus_phi2 = -999;
0348       fclus_theta2 = -999;
0349       fclus_pt2 = -999;
0350       fclus_px2 = -999;
0351       fclus_py2 = -999;
0352       fclus_pz2 = -999;
0353
0354
0355       if(Verbosity() > 1)
0356     {
0357       cout<<"Found one good cluster"<<endl;
0358       cout<<"energy | phi | eta: "<<fclusenergy<<" | "<<fclus_phi<<" | "<<fclus_eta<<endl;
0359     }
0360
0361       //found a first good cluster, lets look for a second
0362       RawClusterContainer::ConstRange fclus2 = fclusters->getClusters();
0363       RawClusterContainer::ConstIterator fclusiter2;
0364       for (fclusiter2 = fclus2.first; fclusiter2 != fclus2.second; ++fclusiter2)
0365       {
0366         RawCluster *cluster2 = fclusiter2->second;
0367
0368         CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetEVec(*cluster2, vertex);
0369         fclusenergy2 = E_vec_cluster2.mag();
0370     //just got the same photon, skip
0371         if (fclusenergy == fclusenergy2)
0372       {
0373         //reset it
0374         fclusenergy2 = -999;
0375         continue;
0376       }
0377
0378         fclus_eta2 = E_vec_cluster2.pseudoRapidity();
0379         fclus_theta2 = E_vec_cluster2.getTheta();
0380         fclus_pt2 = E_vec_cluster2.perp();
0381         fclus_phi2 = E_vec_cluster2.getPhi();
0382
0383
0384
0385         if (fclusenergy2 < mincluspt)
0386           continue;
0387
0388     //avoid double counting, leading photon should always be highest energy
0389         if (fclusenergy < fclusenergy2)
0390           continue;
0391
0392         if (Verbosity() > 1)
0393         {
0394           cout << "identified reco photon 1 energy | phi | eta:  "
0395                << fclusenergy << "  |  " << fclus_phi << "  |  " << fclus_eta << endl;
0396           cout << "identified reco photon 2 energy | phi | eta:  "
0397                << fclusenergy2 << "  |  " << fclus_phi2 << "  |  " << fclus_eta2 << endl;
0398         }
0399
0400         fclus_px2 = fclus_pt2 * TMath::Cos(fclus_phi2);
0401         fclus_py2 = fclus_pt2 * TMath::Sin(fclus_phi2);
0402         fclus_pz2 = fclusenergy2 * TMath::Cos(fclus_theta2);
0403
0404         TLorentzVector phot1, phot2;
0405         phot1.SetPtEtaPhiE(fclus_pt, fclus_eta, fclus_phi, fclusenergy);
0406         phot2.SetPtEtaPhiE(fclus_pt2, fclus_eta2, fclus_phi2, fclusenergy2);
0407
0408         TLorentzVector pi0;
0409         pi0 = phot1 + phot2;
0410
0411     //get the pi0 invmass
0412         invmass = pi0.M();
0413
0414         if (Verbosity() > 1)
0415           cout << "pi0 reco invmass is " << invmass << endl;
0416       }
0417
0418
0419       if(Verbosity() > 1)
0420     {
0421       cout<<"Final clusters found and writing to tree"<<endl;
0422       cout<<"clus 1 - energy | phi | eta: "<<fclusenergy<<" | "<<fclus_phi
0423           <<" | "<<fclus_eta<<endl;
0424       cout<<"clus 2 - energy | phi | eta: "<<fclusenergy2<<" | "<<fclus_phi2
0425           <<" | "<<fclus_eta2<<endl;
0426
0427     }
0428       fcluster_tree->Fill();
0429     }
0430   }
0431
0432
0433
0434
0435
0436
0437   /***********************************************
0438
0439   GET THE CENTRAL ARM EMCAL CLUSTERS
0440
0441   ************************************************/
0442
0443   RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0444   RawClusterContainer::ConstIterator clusiter;
0445
0446   for (clusiter = begin_end.first; clusiter != begin_end.second; ++clusiter)
0447   {
0448     RawCluster *cluster = clusiter->second;
0449
0450     CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0451     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
0452     clus_energy = E_vec_cluster.mag();
0453     clus_eta = E_vec_cluster.pseudoRapidity();
0454     clus_theta = E_vec_cluster.getTheta();
0455     clus_pt = E_vec_cluster.perp();
0456     clus_phi = E_vec_cluster.getPhi();
0457
0458     if (clus_pt < mincluspt)
0459       continue;
0460
0461     TLorentzVector *clus = new TLorentzVector();
0462     clus->SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clus_energy);
0463
0464     float dumarray[4];
0465     clus->GetXYZT(dumarray);
0466     clus_x = dumarray[0];
0467     clus_y = dumarray[1];
0468     clus_z = dumarray[2];
0469     clus_t = dumarray[3];
0470
0471     clus_px = clus_pt * TMath::Cos(clus_phi);
0472     clus_py = clus_pt * TMath::Sin(clus_phi);
0473     clus_pz = sqrt(clus_energy * clus_energy - clus_px * clus_px - clus_py * clus_py);
0474
0475
0476
0477
0478     //reset second cluster variables
0479     clus_energy2 = -999;
0480     clus_eta2 = -999;
0481     clus_theta2 = -999;
0482     clus_pt2 = -999;
0483     clus_phi2 = -999;
0484     clus_px2 = -999;
0485     clus_py2 = -999;
0486     clus_pz2 = -999;
0487
0488
0489
0490
0491
0492     //found a first good cluster, lets look for a second
0493     RawClusterContainer::ConstRange begin_end2 = clusters->getClusters();
0494     RawClusterContainer::ConstIterator clusiter2;
0495
0496     for (clusiter2 = begin_end2.first; clusiter2 != begin_end2.second; ++clusiter2)
0497       {
0498     RawCluster *cluster2 = clusiter2->second;
0499     CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0500     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster2, vertex);
0501     clus_energy2 = E_vec_cluster.mag();
0502     if(clus_energy == clus_energy2)
0503       {
0504         //reset it
0505         clus_energy2 = -999;
0506         continue;
0507       }
0508     clus_eta2 = E_vec_cluster.pseudoRapidity();
0509     clus_theta2 = E_vec_cluster.getTheta();
0510     clus_pt2 = E_vec_cluster.perp();
0511     clus_phi2 = E_vec_cluster.getPhi();
0512
0513     if(clus_energy == clus_energy2 and clus_eta == clus_eta2 and clus_phi == clus_phi2)
0514       continue;
0515
0516     if(clus_energy2 < mincluspt)
0517       continue;
0518
0519     //avoid double counting
0520     if(clus_energy < clus_energy2)
0521       continue;
0522
0523     if (Verbosity() > 1)
0524         {
0525       cout << "CENTRAL ARM ID: "<<endl;
0526           cout << "identified reco photon 1 energy | phi | eta:  "
0527                << clus_energy << "  |  " << clus_phi << "  |  " << clus_eta << endl;
0528           cout << "identified reco photon 2 energy | phi | eta:  "
0529                << clus_energy2 << "  |  " << clus_phi2 << "  |  " << clus_eta2 << endl;
0530         }
0531
0532     clus_px2 = clus_pt2 * TMath::Cos(clus_phi2);
0533     clus_py2 = clus_pt2 * TMath::Sin(clus_phi2);
0534     clus_pz2 = sqrt(clus_energy2 * clus_energy2 - clus_px2 * clus_px2 - clus_py2 * clus_py2);
0535
0536
0537     TLorentzVector phot1, phot2;
0538         phot1.SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clus_energy);
0539         phot2.SetPtEtaPhiE(clus_pt2, clus_eta2, clus_phi2, clus_energy2);
0540
0541         TLorentzVector pi0;
0542         pi0 = phot1 + phot2;
0543
0544     //get the pi0 invmass
0545         cent_invmass = pi0.M();
0546
0547         if (Verbosity() > 1)
0548           cout << "Central pi0 reco invmass is " << cent_invmass << endl;
0549
0550
0551     cluster_tree->Fill();
0552       }
0553
0554   }
0555
0556   tree->Fill();
0557
0558
0559
0560   return 0;
0561 }
0562
0563 int Forward_pi0s::End(PHCompositeNode *topnode)
0564 {
0565   std::cout << " DONE PROCESSING " << endl;
0566
0567   file->Write();
0568   file->Close();
0569   return 0;
0570 }
0571
0572 void Forward_pi0s::Set_Tree_Branches()
0573 {
0574   inhcal_tree = new TTree("hcalclustree", "a tree with inner hcal cluster information");
0575   inhcal_tree->Branch("hcal_energy", &hcal_energy, "hcal_energy/F");
0576   inhcal_tree->Branch("hcal_eta", &hcal_eta, "hcal_eta/F");
0577   inhcal_tree->Branch("hcal_phi", &hcal_phi, "hcal_phi/F");
0578   inhcal_tree->Branch("hcal_pt", &hcal_pt, "hcal_pt/F");
0579   inhcal_tree->Branch("hcal_theta", &hcal_theta, "hcal_theta/F");
0580   inhcal_tree->Branch("hcal_x", &hcal_x, "hcal_x/F");
0581   inhcal_tree->Branch("hcal_y", &hcal_y, "hcal_y/F");
0582   inhcal_tree->Branch("hcal_z", &hcal_z, "hcal_z/F");
0583   inhcal_tree->Branch("hcal_t", &hcal_t, "hcal_t/F");
0584   inhcal_tree->Branch("hcal_px", &hcal_px, "hcal_px/F");
0585   inhcal_tree->Branch("hcal_py", &hcal_py, "hcal_py/F");
0586   inhcal_tree->Branch("hcal_pz", &hcal_pz, "hcal_pz/F");
0587   inhcal_tree->Branch("nevents", &nevents, "nevents/I");
0588   inhcal_tree->Branch("hclustruthpx", &hclustruthpx, "hclustruthpx/F");
0589   inhcal_tree->Branch("hclustruthpy", &hclustruthpy, "hclustruthpy/F");
0590   inhcal_tree->Branch("hclustruthpz", &hclustruthpz, "hclustruthpz/F");
0591   inhcal_tree->Branch("hclustruthenergy", &hclustruthenergy, "hclustruthenergy/F");
0592   inhcal_tree->Branch("hclustruthpt", &hclustruthpt, "hclustruthpt/F");
0593   inhcal_tree->Branch("hclustruthphi", &hclustruthphi, "hclustruthphi/F");
0594   inhcal_tree->Branch("hclustrutheta", &hclustrutheta, "hclustrutheta/F");
0595
0596   cluster_tree = new TTree("clustertree", "A tree with EMCal cluster information");
0597   cluster_tree->Branch("clus_energy", &clus_energy, "clus_energy/F");
0598   cluster_tree->Branch("clus_eta", &clus_eta, "clus_eta/F");
0599   cluster_tree->Branch("clus_phi", &clus_phi, "clus_phi/F");
0600   cluster_tree->Branch("clus_pt", &clus_pt, "clus_pt/F");
0601   cluster_tree->Branch("clus_theta", &clus_theta, "clus_theta/F");
0602   cluster_tree->Branch("clus_x", &clus_x, "clus_x/F");
0603   cluster_tree->Branch("clus_y", &clus_y, "clus_y/F");
0604   cluster_tree->Branch("clus_z", &clus_z, "clus_z/F");
0605   cluster_tree->Branch("clus_t", &clus_t, "clus_t/F");
0606   cluster_tree->Branch("clus_px", &clus_px, "clus_px/F");
0607   cluster_tree->Branch("clus_py", &clus_py, "clus_py/F");
0608   cluster_tree->Branch("clus_pz", &clus_pz, "clus_pz/F");
0609   cluster_tree->Branch("nevents", &nevents, "nevents/I");
0610
0611   cluster_tree->Branch("clus_px2", &clus_px2, "clus_px2/F");
0612   cluster_tree->Branch("clus_py2", &clus_py2, "clus_py2/F");
0613   cluster_tree->Branch("clus_pz2", &clus_pz2, "clus_pz2/F");
0614   cluster_tree->Branch("clus_energy2", &clus_energy2, "clus_energy2/F");
0615   cluster_tree->Branch("clus_eta2", &clus_eta2, "clus_eta2/F");
0616   cluster_tree->Branch("clus_phi2", &clus_phi2, "clus_phi2/F");
0617   cluster_tree->Branch("clus_pt2", &clus_pt2, "clus_pt2/F");
0618   cluster_tree->Branch("clus_theta2", &clus_theta2, "clus_theta2/F");
0619   cluster_tree->Branch("Cent_invmass", &cent_invmass,"cent_invmass/F");
0620
0621   cluster_tree->Branch("tpi0e", &tpi0e, "tpi0e/F");
0622   cluster_tree->Branch("tpi0px", &tpi0px, "tpi0px/F");
0623   cluster_tree->Branch("tpi0py", &tpi0py, "tpi0py/F");
0624   cluster_tree->Branch("tpi0pz", &tpi0pz, "tpi0pz/F");
0625   cluster_tree->Branch("tpi0pid", &tpi0pid, "tpi0pid/F");
0626   cluster_tree->Branch("tpi0pt", &tpi0pt, "tpi0pt/F");
0627   cluster_tree->Branch("tpi0phi", &tpi0phi, "tpi0phi/F");
0628   cluster_tree->Branch("tpi0eta", &tpi0eta, "tpi0eta/F");
0629   cluster_tree->Branch("tphote1", &tphote1, "tphote1/F");
0630   cluster_tree->Branch("tphotpx1", &tphotpx1, "tphotpx1/F");
0631   cluster_tree->Branch("tphotpy1", &tphotpy1, "tphotpy1/F");
0632   cluster_tree->Branch("tphotpz1", &tphotpz1, "tphotpz1/F");
0633   cluster_tree->Branch("tphotpt1", &tphotpt1, "tphotpt1/F");
0634   cluster_tree->Branch("tphotpid1", &tphotpid1, "tphotpid1/I");
0635   cluster_tree->Branch("tphotparentid1", &tphotparentid1, "tphotparentid1/I");
0636   cluster_tree->Branch("tphotphi1", &tphotphi1, "tphotphi1/F");
0637   cluster_tree->Branch("tphoteta1", &tphoteta1, "tphoteta1/F");
0638
0639   cluster_tree->Branch("tphote2", &tphote2, "tphote2/F");
0640   cluster_tree->Branch("tphotpx2", &tphotpx2, "tphotpx2/F");
0641   cluster_tree->Branch("tphotpy2", &tphotpy2, "tphotpy2/F");
0642   cluster_tree->Branch("tphotpz2", &tphotpz2, "tphotpz2/F");
0643   cluster_tree->Branch("tphotpt2", &tphotpt2, "tphotpt2/F");
0644   cluster_tree->Branch("tphotpid2", &tphotpid2, "tphotpid2/I");
0645   cluster_tree->Branch("tphotparentid2", &tphotparentid2, "tphotparentid2/I");
0646   cluster_tree->Branch("tphotphi2", &tphotphi2, "tphotphi2/F");
0647   cluster_tree->Branch("tphoteta2", &tphoteta2, "tphoteta2/F");
0648
0649
0650
0651   fcluster_tree = new TTree("fclustertree", "A tree with FEMCal cluster information");
0652   fcluster_tree->Branch("invmass", &invmass, "invmass/F");
0653   fcluster_tree->Branch("fclusenergy", &fclusenergy, "fclusenergy/F");
0654   fcluster_tree->Branch("fclus_eta", &fclus_eta, "fclus_eta/F");
0655   fcluster_tree->Branch("fclus_phi", &fclus_phi, "fclus_phi/F");
0656   fcluster_tree->Branch("fclus_pt", &fclus_pt, "fclus_pt/F");
0657   fcluster_tree->Branch("fclus_theta", &fclus_theta, "fclus_theta/F");
0658   fcluster_tree->Branch("fclus_px", &fclus_px, "fclus_px/F");
0659   fcluster_tree->Branch("fclus_py", &fclus_py, "fclus_py/F");
0660   fcluster_tree->Branch("fclus_pz", &fclus_pz, "fclus_pz/F");
0661   fcluster_tree->Branch("nevents", &nevents, "nevents/I");
0662   fcluster_tree->Branch("fclusenergy2", &fclusenergy2, "fclusenergy2/F");
0663   fcluster_tree->Branch("fclus_eta2", &fclus_eta2, "fclus_eta2/F");
0664   fcluster_tree->Branch("fclus_phi2", &fclus_phi2, "fclus_phi2/F");
0665   fcluster_tree->Branch("fclus_pt2", &fclus_pt2, "fclus_pt2/F");
0666   fcluster_tree->Branch("fclus_theta2", &fclus_theta2, "fclus_theta2/F");
0667   fcluster_tree->Branch("fclus_px2", &fclus_px2, "fclus_px2/F");
0668   fcluster_tree->Branch("fclus_py2", &fclus_py2, "fclus_py2/F");
0669   fcluster_tree->Branch("fclus_pz2", &fclus_pz2, "fclus_pz2/F");
0670
0671   fcluster_tree->Branch("tpi0e", &tpi0e, "tpi0e/F");
0672   fcluster_tree->Branch("tpi0px", &tpi0px, "tpi0px/F");
0673   fcluster_tree->Branch("tpi0py", &tpi0py, "tpi0py/F");
0674   fcluster_tree->Branch("tpi0pz", &tpi0pz, "tpi0pz/F");
0675   fcluster_tree->Branch("tpi0pid", &tpi0pid, "tpi0pid/F");
0676   fcluster_tree->Branch("tpi0pt", &tpi0pt, "tpi0pt/F");
0677   fcluster_tree->Branch("tpi0phi", &tpi0phi, "tpi0phi/F");
0678   fcluster_tree->Branch("tpi0eta", &tpi0eta, "tpi0eta/F");
0679   fcluster_tree->Branch("tphote1", &tphote1, "tphote1/F");
0680   fcluster_tree->Branch("tphotpx1", &tphotpx1, "tphotpx1/F");
0681   fcluster_tree->Branch("tphotpy1", &tphotpy1, "tphotpy1/F");
0682   fcluster_tree->Branch("tphotpz1", &tphotpz1, "tphotpz1/F");
0683   fcluster_tree->Branch("tphotpt1", &tphotpt1, "tphotpt1/F");
0684   fcluster_tree->Branch("tphotpid1", &tphotpid1, "tphotpid1/I");
0685   fcluster_tree->Branch("tphotparentid1", &tphotparentid1, "tphotparentid1/I");
0686   fcluster_tree->Branch("tphotphi1", &tphotphi1, "tphotphi1/F");
0687   fcluster_tree->Branch("tphoteta1", &tphoteta1, "tphoteta1/F");
0688
0689   fcluster_tree->Branch("tphote2", &tphote2, "tphote2/F");
0690   fcluster_tree->Branch("tphotpx2", &tphotpx2, "tphotpx2/F");
0691   fcluster_tree->Branch("tphotpy2", &tphotpy2, "tphotpy2/F");
0692   fcluster_tree->Branch("tphotpz2", &tphotpz2, "tphotpz2/F");
0693   fcluster_tree->Branch("tphotpt2", &tphotpt2, "tphotpt2/F");
0694   fcluster_tree->Branch("tphotpid2", &tphotpid2, "tphotpid2/I");
0695   fcluster_tree->Branch("tphotparentid2", &tphotparentid2, "tphotparentid2/I");
0696   fcluster_tree->Branch("tphotphi2", &tphotphi2, "tphotphi2/F");
0697   fcluster_tree->Branch("tphoteta2", &tphoteta2, "tphoteta2/F");
0698
0699   truth_g4particles = new TTree("truthtree_g4", "a tree with all truth g4 particles");
0700   truth_g4particles->Branch("truthpx", &truthpx, "truthpx/F");
0701   truth_g4particles->Branch("truthpy", &truthpy, "truthpy/F");
0702   truth_g4particles->Branch("truthpz", &truthpz, "truthpz/F");
0703   truth_g4particles->Branch("truthp", &truthp, "truthp/F");
0704   truth_g4particles->Branch("truthenergy", &truthenergy, "truthenergy/F");
0705   truth_g4particles->Branch("truthphi", &truthphi, "truthphi/F");
0706   truth_g4particles->Branch("trutheta", &trutheta, "trutheta/F");
0707   truth_g4particles->Branch("truthpt", &truthpt, "truthpt/F");
0708   truth_g4particles->Branch("truthpid", &truthpid, "truthpid/I");
0709   truth_g4particles->Branch("nevents", &nevents, "nevents/I");
0710
0711   truth_g4particles->Branch("tphote1", &tphote1, "tphote1/F");
0712   truth_g4particles->Branch("tphotpx1", &tphotpx1, "tphotpx1/F");
0713   truth_g4particles->Branch("tphotpy1", &tphotpy1, "tphotpy1/F");
0714   truth_g4particles->Branch("tphotpz1", &tphotpz1, "tphotpz1/F");
0715   truth_g4particles->Branch("tphotpt1", &tphotpt1, "tphotpt1/F");
0716   truth_g4particles->Branch("tphotpid1", &tphotpid1, "tphotpid1/I");
0717   truth_g4particles->Branch("tphotparentid1", &tphotparentid1, "tphotparentid1/I");
0718   truth_g4particles->Branch("tphotphi1", &tphotphi1, "tphotphi1/F");
0719   truth_g4particles->Branch("tphoteta1", &tphoteta1, "tphoteta1/F");
0720
0721   truth_g4particles->Branch("tphote2", &tphote2, "tphote2/F");
0722   truth_g4particles->Branch("tphotpx2", &tphotpx2, "tphotpx2/F");
0723   truth_g4particles->Branch("tphotpy2", &tphotpy2, "tphotpy2/F");
0724   truth_g4particles->Branch("tphotpz2", &tphotpz2, "tphotpz2/F");
0725   truth_g4particles->Branch("tphotpt2", &tphotpt2, "tphotpt2/F");
0726   truth_g4particles->Branch("tphotpid2", &tphotpid2, "tphotpid2/I");
0727   truth_g4particles->Branch("tphotparentid2", &tphotparentid2, "tphotparentid2/I");
0728   truth_g4particles->Branch("tphotphi2", &tphotphi2, "tphotphi2/F");
0729   truth_g4particles->Branch("tphoteta2", &tphoteta2, "tphoteta2/F");
0730 }
0731
0732 void Forward_pi0s::initialize_things()
0733 {
0734   hcal_energy = -999;
0735   hcal_eta = -999;
0736   hcal_phi = -999;
0737   hcal_pt = -999;
0738   hcal_px = -999;
0739   hcal_py = -999;
0740   hcal_pz = -999;
0741   hcal_theta = -999;
0742   hcal_x = -999;
0743   hcal_y = -999;
0744   hcal_z = -999;
0745   hcal_t = -999;
0746
0747   clus_energy = -999;
0748   clus_eta = -999;
0749   clus_phi = -999;
0750   clus_pt = -999;
0751   clus_px = -999;
0752   clus_py = -999;
0753   clus_pz = -999;
0754   clus_theta = -999;
0755   clus_x = -999;
0756   clus_y = -999;
0757   clus_z = -999;
0758   clus_t = -999;
0759
0760   fclusenergy2 = -999;
0761   fclus_eta2 = -999;
0762   fclus_phi2 = -999;
0763   fclus_theta2 = -999;
0764   fclus_pt2 = -999;
0765   fclus_px2 = -999;
0766   fclus_py2 = -999;
0767   fclus_pz2 = -999;
0768   invmass = -999;
0769
0770   tphote1 = -999;
0771   tphotpx1 = -999;
0772   tphotpy1 = -999;
0773   tphotpz1 = -999;
0774   tphotpid1 = -999;
0775   tphotparentid1 = -999;
0776   tphotpt1 = -999;
0777   tphotphi1 = -999;
0778   tphoteta1 = -999;
0779   tphote2 = -999;
0780   tphotpx2 = -999;
0781   tphotpy2 = -999;
0782   tphotpz2 = -999;
0783   tphotpid2 = -999;
0784   tphotparentid2 = -999;
0785   tphotpt2 = -999;
0786   tphotphi2 = -999;
0787   tphoteta2 = -999;
0788   fclusenergy = -999;
0789   fclus_eta = -999;
0790   fclus_phi = -999;
0791   fclus_theta = -999;
0792   fclus_pt = -999;
0793   fclustruthpid = -999;
0794   fclustruthpx = -999;
0795   fclustruthpy = -999;
0796   fclustruthpz = -999;
0797   fclustruthenergy = -999;
0798   fclustruthpt = -999;
0799   fclustruthphi = -999;
0800   fclustrutheta = -999;
0801   fclus_px = -999;
0802   fclus_py = -999;
0803   fclus_pz = -999;
0804
0805   tpi0e = -999;
0806   tpi0px = -999;
0807   tpi0py = -999;
0808   tpi0pz = -999;
0809   tpi0pid = -999;
0810   tpi0pt = -999;
0811   tpi0phi = -999;
0812   tpi0eta = -999;
0813 }