EventMix/ana_calo/sPHAnalysis_calo.cc

0001 #include "sPHAnalysis_calo.h"
0002
0003 #include <cmath>
0004 #include <cstdlib>
0005 #include <cstdio>
0006 #include <iostream>
0007 #include <iomanip>
0008 #include <fstream>
0009
0010 #include "TFile.h"
0011 #include "TNtuple.h"
0012 #include "TH1F.h"
0013 #include "TH1D.h"
0014 #include "TF1.h"
0015 #include "TLorentzVector.h"
0016 #include "TRandom2.h"
0017
0018 #include <globalvertex/GlobalVertexMap.h>
0019 #include <globalvertex/GlobalVertex.h>
0020
0021 #include <calobase/RawClusterContainer.h>
0022 #include <calobase/RawCluster.h>
0023 #include <calobase/RawClusterv1.h>
0024 #include <calobase/RawTowerv2.h>
0025 #include <calobase/RawTowerGeom.h>
0026 #include <calobase/RawTowerContainer.h>
0027 #include <calobase/TowerInfov1.h>
0028 #include <calobase/TowerInfoContainerv1.h>
0029 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0030 #include <calobase/RawTowerGeomContainer.h>
0031
0032 #include <phool/getClass.h>
0033 #include <phool/recoConsts.h>
0034 #include <phool/PHCompositeNode.h>
0035 #include <phool/PHIODataNode.h>
0036 #include <phool/PHRandomSeed.h>
0037 #include <fun4all/Fun4AllReturnCodes.h>
0038
0039 #include <gsl/gsl_rng.h>
0040
0041 using namespace std;
0042
0043 //==============================================================
0044
0045 sPHAnalysis_calo::sPHAnalysis_calo(const std::string &name, const std::string &filename) : SubsysReco(name)
0046 {
0047   OutputNtupleFile=nullptr;
0048   OutputFileName=filename;
0049   EventNumber=0;
0050   ntp_notracking=nullptr;
0051   h_mult=nullptr;
0052   h_ecore= nullptr;
0053   _whattodo = 0;
0054   _rng = nullptr;
0055 }
0056
0057 //==============================================================
0058
0059 int sPHAnalysis_calo::Init(PHCompositeNode *topNode)
0060 {
0061   std::cout << "sPHAnalysis_calo::Init started..." << endl;
0062   OutputNtupleFile = new TFile(OutputFileName.c_str(),"RECREATE");
0063   std::cout << "sPHAnalysis_calo::Init: output file " << OutputFileName.c_str() << " opened." << endl;
0064
0065   ntp_notracking = new  TNtuple("ntp_notracking","","mass:pt:e1:e2:hoe1:hoe2");
0066   h_mult =  new TH1D("h_mult","",1000,0.,10000.);
0067   h_ecore =  new TH1D("h_ecore","",200,0.,20.);
0068
0069   char hname[99];
0070   for(int i=0; i<256*96; i++) {
0071     sprintf(hname,"h_pedestal_%d",i);
0072     h_pedestal[i] = new TH1F(hname,hname,10000,-0.5,9999.5);
0073   }
0074
0075   _rng = new TRandom2();
0076   _rng->SetSeed(0);
0077
0078   cout << "sPHAnalysis_calo::Init ended." << endl;
0079   return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 //==============================================================
0083
0084 int sPHAnalysis_calo::InitRun(PHCompositeNode *topNode)
0085 {
0086   cout << "sPHAnalysis_calo::InitRun started..." << endl;
0087   cout << "sPHAnalysis_calo::InitRun ended." << endl;
0088   return Fun4AllReturnCodes::EVENT_OK;
0089 }
0090
0091 //==============================================================
0092
0093 int sPHAnalysis_calo::process_event(PHCompositeNode *topNode)
0094 {
0095   if(_whattodo==0) {
0096     return process_event_test(topNode); // simulation
0097   } else if(_whattodo==1) {
0098     return process_event_data(topNode); // data
0099   } else {
0100     cerr << "ERROR: wrong choice of what to do." << endl; return Fun4AllReturnCodes::ABORTRUN;
0101   }
0102 }
0103
0104 //======================================================================
0105
0106 int sPHAnalysis_calo::process_event_data(PHCompositeNode *topNode) {
0107   EventNumber++;
0108   float tmp1[99];
0109
0110   cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
0111   if(EventNumber==1) topNode->print();
0112
0113   TowerInfoContainer* _towersCEMC = findNode::getClass<TowerInfoContainerv1>(topNode, "TOWERS_CEMC");
0114   if (!_towersCEMC) { std::cerr<<"No TOWERS_CEMC Node"<<std::endl; return Fun4AllReturnCodes::ABORTEVENT; }
0115   unsigned int ntowers = _towersCEMC->size();
0116   cout << " # of towers = " << ntowers << endl;
0117
0118   for (unsigned int channel = 0; channel < ntowers; channel++)
0119   {
0120     TowerInfo *cemcinfo_raw = _towersCEMC->get_tower_at_channel(channel);
0121     float raw_amplitude = cemcinfo_raw->get_energy();
0122     float raw_time = cemcinfo_raw->get_time();
0123     unsigned int key = _towersCEMC->encode_key(channel);
0124     unsigned int iphi = _towersCEMC->getTowerPhiBin(key);
0125     unsigned int ieta = _towersCEMC->getTowerEtaBin(key);
0126     h_pedestal[channel]->Fill(raw_amplitude);
0127     if(channel>0 && channel<20) cout << channel << " " << iphi << " " << ieta << " " << raw_amplitude << " " << raw_time << endl;
0128   }
0129
0130   cout << "sPHAnalysis_calo::procedss_event_data() ended." << endl;
0131   return Fun4AllReturnCodes::EVENT_OK;
0132 }
0133
0134 //======================================================================
0135
0136 int sPHAnalysis_calo::process_event_test(PHCompositeNode *topNode) {
0137   EventNumber++;
0138   float tmp1[99];
0139
0140   cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
0141   if(EventNumber==1) topNode->print();
0142
0143 // Find event vertex
0144   GlobalVertexMap *global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0145   //cout << "Number of GlobalVertexMap entries = " << global_vtxmap->size() << endl;
0146   double Zvtx = 0.;
0147   for (GlobalVertexMap::Iter iter = global_vtxmap->begin(); iter != global_vtxmap->end(); ++iter)
0148   {
0149     GlobalVertex *vtx = iter->second;
0150     if(vtx->get_id()==1) Zvtx = vtx->get_z(); // BBC vertex (?)
0151     //cout << "Global vertex: " << vtx->get_id() << " " << vtx->get_z() << " " << vtx->get_t() << endl;
0152   }
0153   cout << "Global vertex Z = " << Zvtx << endl;
0154
0155 // Find node with CEMC clusters
0156   RawClusterContainer* cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0157   if(!cemc_clusters) {
0158     cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
0159     return Fun4AllReturnCodes::ABORTEVENT;
0160   }
0161   cout << "   Number of CEMC clusters = " << cemc_clusters->size() << endl;
0162   h_mult->Fill(cemc_clusters->size());
0163
0164 // Fetch CEMC geometry
0165   RawTowerGeomContainer* _geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0166   if(!_geomEM) std::cerr<<"No TOWERGEOM_CEMC"<<std::endl;
0167   RawTowerGeomContainer* _geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0168   if(!_geomIH) std::cerr<<"No TOWERGEOM_HCALIN"<<std::endl;
0169   RawTowerGeomContainer* _geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0170   if(!_geomIH) std::cerr<<"No TOWERGEOM_HCALIN"<<std::endl;
0171
0172   RawTowerContainer* _towersRawEM = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0173   if (!_towersRawEM) std::cerr<<"No TOWER_CALIB_CEMC Node"<<std::endl;
0174   RawTowerContainer* _towersRawIH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0175   if (!_towersRawIH) std::cerr<<"No TOWER_CALIB_HCALIN Node"<<std::endl;
0176   RawTowerContainer* _towersRawOH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0177   if (!_towersRawOH) std::cerr<<"No TOWER_CALIB_HCALOUT Node"<<std::endl;
0178
0179   vector<TLorentzVector> electrons;
0180   vector<double> vhoe;
0181
0182     // loop over CEMC clusters with E > 2 GeV
0183     RawClusterContainer::Range begin_end = cemc_clusters->getClusters();
0184     RawClusterContainer::Iterator iter;
0185     for (iter = begin_end.first; iter != begin_end.second; ++iter)
0186     {
0187       RawCluster* cluster = iter->second;
0188       if(!cluster) { cout << "ERROR: bad cluster pointer = " << cluster << endl; continue; }
0189       else {
0190         double cemc_ecore = cluster->get_ecore();
0191         h_ecore->Fill(cemc_ecore);
0192         if(cemc_ecore<2.0) continue;
0193         double cemc_x = cluster->get_x();
0194         double cemc_y = cluster->get_y();
0195         double cemc_z = cluster->get_z() - Zvtx; // correct for event vertex position
0196         //double cemc_r = cluster->get_r();
0197         double cemc_r = sqrt(pow(cemc_x,2)+pow(cemc_y,2));
0198         double cemc_rr = sqrt(pow(cemc_r,2)+pow(cemc_z,2));
0199         double cemc_eta = asinh( cemc_z / cemc_r );
0200         double cemc_phi = atan2( cemc_y, cemc_x );
0201         double cemc_px = cemc_ecore * (cemc_x/cemc_rr);
0202         double cemc_py = cemc_ecore * (cemc_y/cemc_rr);
0203         double cemc_pz = cemc_ecore * (cemc_z/cemc_rr);
0204         //cout << "CEMC cluster: " << cemc_ecore << " " << cemc_eta << " " << cemc_phi << endl;
0205
0206         // find closest CEMC tower in eta-phi space
0207           double distem = 99999.;
0208           RawTower* thetowerem = nullptr;
0209           RawTowerContainer::ConstRange begin_end_rawEM = _towersRawEM->getTowers();
0210           for (RawTowerContainer::ConstIterator rtiter = begin_end_rawEM.first; rtiter != begin_end_rawEM.second; ++rtiter) {
0211             RawTower *tower = rtiter->second;
0212             RawTowerGeom *tower_geom = _geomEM->get_tower_geometry(tower->get_key());
0213             double cemc_tower_phi  = tower_geom->get_phi();
0214             double cemc_tower_x = tower_geom->get_center_x();
0215             double cemc_tower_y = tower_geom->get_center_y();
0216             double cemc_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
0217             double cemc_tower_r = sqrt(pow(cemc_tower_x,2)+pow(cemc_tower_y,2));
0218             double cemc_tower_eta = asinh( cemc_tower_z / cemc_tower_r );
0219             double tmpdist = sqrt(pow(cemc_eta-cemc_tower_eta,2)+pow(cemc_phi-cemc_tower_phi,2));
0220             if(tmpdist<distem) { distem = tmpdist; thetowerem = tower; }
0221           }
0222           RawTowerGeom *thetower_geom_em = _geomEM->get_tower_geometry(thetowerem->get_key());
0223           unsigned int ietaem = thetower_geom_em->get_bineta();
0224           unsigned int jphiem = thetower_geom_em->get_binphi();
0225
0226           // find closest HCALIN tower in eta-phi space
0227           double distin = 99999.;
0228           RawTower* thetowerin = nullptr;
0229           RawTowerContainer::ConstRange begin_end_rawIN = _towersRawIH->getTowers();
0230           for (RawTowerContainer::ConstIterator rtiter = begin_end_rawIN.first; rtiter != begin_end_rawIN.second; ++rtiter) {
0231             RawTower *tower = rtiter->second;
0232             RawTowerGeom *tower_geom = _geomIH->get_tower_geometry(tower->get_key());
0233             double hcalin_tower_phi  = tower_geom->get_phi();
0234             double hcalin_tower_x = tower_geom->get_center_x();
0235             double hcalin_tower_y = tower_geom->get_center_y();
0236             double hcalin_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
0237             double hcalin_tower_r = sqrt(pow(hcalin_tower_x,2)+pow(hcalin_tower_y,2));
0238             double hcalin_tower_eta = asinh( hcalin_tower_z / hcalin_tower_r );
0239             double tmpdist = sqrt(pow(cemc_eta-hcalin_tower_eta,2)+pow(cemc_phi-hcalin_tower_phi,2));
0240             if(tmpdist<distin) { distin = tmpdist; thetowerin = tower; }
0241           }
0242           RawTowerGeom *thetower_geom = _geomIH->get_tower_geometry(thetowerin->get_key());
0243           unsigned int ieta = thetower_geom->get_bineta();
0244           unsigned int jphi = thetower_geom->get_binphi();
0245
0246           // Calcuate 3x3 energy deposit in HCALIN. Inner HCAL has 24 bins in eta and 64 bins in phi
0247           double e3x3in = 0.;
0248           if(ieta<1 || ieta>22) continue; // ignore clusters in edge towers
0249           for(unsigned int i=0; i<=2; i++) {
0250             for(unsigned int j=0; j<=2; j++) {
0251               unsigned int itmp = ieta-1+i;
0252               unsigned int jtmp = 0;
0253               if(jphi==0 && j==0) { jtmp = 63; } // wrap around
0254               else if(jphi==63 && j==2) { jtmp = 0; } // wrap around
0255               else { jtmp = jphi-1+j; }
0256               RawTower* tmptower = _towersRawIH->getTower(itmp,jtmp);
0257               if(tmptower) { e3x3in += tmptower->get_energy(); }
0258             }
0259           }
0260           //h_notracking_eoh->Fill(e3x3in/cemc_ecore);
0261
0262           // find closest HCALOUT tower in eta-phi space
0263           double distout = 99999.;
0264           RawTower* thetowerout = nullptr;
0265           RawTowerContainer::ConstRange begin_end_rawON = _towersRawOH->getTowers();
0266           for (RawTowerContainer::ConstIterator rtiter = begin_end_rawON.first; rtiter != begin_end_rawON.second; ++rtiter) {
0267             RawTower *tower = rtiter->second;
0268             RawTowerGeom *tower_geom = _geomOH->get_tower_geometry(tower->get_key());
0269             double hcalout_tower_phi  = tower_geom->get_phi();
0270             double hcalout_tower_x = tower_geom->get_center_x();
0271             double hcalout_tower_y = tower_geom->get_center_y();
0272             double hcalout_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
0273             double hcalout_tower_r = sqrt(pow(hcalout_tower_x,2)+pow(hcalout_tower_y,2));
0274             double hcalout_tower_eta = asinh( hcalout_tower_z / hcalout_tower_r );
0275             double tmpdist = sqrt(pow(cemc_eta-hcalout_tower_eta,2)+pow(cemc_phi-hcalout_tower_phi,2));
0276             if(tmpdist<distout) { distout = tmpdist; thetowerout = tower; }
0277           }
0278           RawTowerGeom *thetower_geomout = _geomOH->get_tower_geometry(thetowerout->get_key());
0279           unsigned int ietaout = thetower_geomout->get_bineta();
0280           unsigned int jphiout = thetower_geomout->get_binphi();
0281
0282           // Calcuate 3x3 energy deposit in HCALOUT. Outer HCAL has 24 bins in eta and 64 bins in phi
0283           double e3x3out = 0.;
0284           if(ietaout<1 || ietaout>22) continue; // ignore clusters in edge towers
0285           for(unsigned int i=0; i<=2; i++) {
0286             for(unsigned int j=0; j<=2; j++) {
0287               unsigned int itmp = ietaout-1+i;
0288               unsigned int jtmp = 0;
0289               if(jphiout==0 && j==0) { jtmp = 63; } // wrap around
0290               else if(jphiout==63 && j==2) { jtmp = 0; } // wrap around
0291               else { jtmp = jphiout-1+j; }
0292               RawTower* tmptower = _towersRawOH->getTower(itmp,jtmp);
0293               if(tmptower) { e3x3out += tmptower->get_energy(); }
0294             }
0295           }
0296
0297           if(e3x3in/cemc_ecore>0.06) continue; // reject hadrons, 90% eID efficiency is with 0.058 cut
0298
0299           TLorentzVector tmp = TLorentzVector(cemc_px,cemc_py,cemc_pz,cemc_ecore);
0300           electrons.push_back(tmp);
0301           vhoe.push_back(e3x3in/cemc_ecore);
0302
0303       } // valid CEMC cluster
0304     } // end loop over CEMC clusters
0305     cout << "number of selected electrons = " << electrons.size() << " " << vhoe.size() << endl;
0306
0307 // Make Upsilons
0308
0309   if(electrons.size()>=2) {
0310     for(long unsigned int i=0; i<electrons.size()-1; i++) {
0311       for(long unsigned int j=i+1; j<electrons.size(); j++) {
0312         TLorentzVector pair = electrons[i]+electrons[j];
0313         cout << i << " " << j << endl;
0314         tmp1[0] = pair.M();
0315         tmp1[1] = pair.Pt();
0316         cout << pair.M() << " " << pair.Pt() << endl;
0317         tmp1[2] = (electrons[i]).Pt();
0318         tmp1[3] = (electrons[j]).Pt();
0319         cout << vhoe[i] << " " << vhoe[j] << endl;
0320         tmp1[4] = vhoe[i];
0321         tmp1[5] = vhoe[j];
0322         cout << "filling..." << endl;
0323           ntp_notracking->Fill(tmp1);
0324         cout << "done." << endl;
0325       }
0326     }
0327   }
0328
0329   cout << "sPHAnalysis_calo::procedss_event_test() ended." << endl;
0330   return Fun4AllReturnCodes::EVENT_OK;
0331 }
0332
0333 //======================================================================
0334
0335 int sPHAnalysis_calo::End(PHCompositeNode *topNode)
0336 {
0337   std::cout << "sPHAnalysis_calo::End() started, Number of processed events = " << EventNumber << endl;
0338   OutputNtupleFile->Write();
0339   OutputNtupleFile->Close();
0340   return Fun4AllReturnCodes::EVENT_OK;
0341 }
0342
0343