Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:06 

0001 #include "ShowerSize.h"
0002 
0003 #include <g4main/PHG4HitContainer.h>
0004 #include <g4main/PHG4Hit.h>
0005 
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4VtxPoint.h>
0009 
0010 #include <fun4all/Fun4AllHistoManager.h>
0011 
0012 #include <phool/getClass.h>
0013 
0014 #include <TFile.h>
0015 #include <TH1.h>
0016 #include <TH2.h>
0017 #include <TNtuple.h>
0018 
0019 #include<sstream>
0020 
0021 using namespace std;
0022 
0023 ShowerSize::ShowerSize(const std::string &name, const std::string &filename):
0024   SubsysReco( name ),
0025   nblocks(0),
0026   hm(nullptr),
0027   _filename(filename),
0028   ntups(nullptr),
0029   ntupe(nullptr),
0030   ntup(nullptr),
0031   outfile(nullptr)
0032 {}
0033 
0034 ShowerSize::~ShowerSize()
0035 {
0036   //  delete ntup;
0037   delete hm;
0038 } 
0039 
0040 
0041 int
0042 ShowerSize::Init( PHCompositeNode* )
0043 {
0044   ostringstream hname, htit;
0045   hm = new  Fun4AllHistoManager(Name());
0046   outfile = new TFile(_filename.c_str(), "RECREATE");
0047   ntups = new TNtuple("sz"," Shower size","rad:em:hi:ho:mag:bh");
0048   ntupe = new TNtuple("truth", "The Absolute Truth", "phi:theta:eta:e:p");
0049   ntup = new TNtuple("de", "Change in Angles", "ID:dphi:dtheta:dtotal:edep");
0050   return 0;
0051 }
0052 
0053 int
0054 ShowerSize::process_event( PHCompositeNode* topNode )
0055 {
0056   map<int, PHG4Particle*>::const_iterator particle_iter;
0057   PHG4TruthInfoContainer *_truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0058 
0059   PHG4TruthInfoContainer::ConstRange primary_range =
0060     _truth_container->GetPrimaryParticleRange();
0061   float ntvars[5] = {0};
0062   for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
0063        particle_iter != primary_range.second; ++particle_iter)
0064   {
0065     PHG4Particle *particle = particle_iter->second;
0066     ntvars[0] = atan2(particle->get_py(), particle->get_px());
0067     ntvars[1] = atan(sqrt(particle->get_py()*particle->get_py() + 
0068               particle->get_px()*particle->get_px()) /
0069              particle->get_pz());
0070     if (ntvars[1] < 0)
0071       {
0072     ntvars[1]+=M_PI;
0073       }
0074     ntvars[2] = 0.5*log((particle->get_e()+particle->get_pz())/
0075             (particle->get_e()-particle->get_pz()));
0076     ntvars[3] = particle->get_e();
0077     ntvars[4] = particle->get_pz();
0078   }
0079   ntupe->Fill(ntvars);
0080   PHG4HitContainer *hits = nullptr;
0081   double emc[10] = {0};
0082   double ih[10] = {0};
0083   double oh[10] = {0};
0084   double mag[10] = {0};
0085   double bh[10] = {0};
0086   double eall[5] = {0};
0087   float detid[5] = {0};
0088   string nodename[2] = {"G4HIT_CEMC","G4HIT_ABSORBER_CEMC"};
0089   for (int j=0; j<2;j++)
0090   {
0091     hits = findNode::getClass<PHG4HitContainer>(topNode,nodename[j]);
0092     if (hits)
0093     {
0094       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0095       for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0096       {
0097     // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0098     //           hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0099     //           hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0100     double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0101     double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0102                  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0103                 hit_iter->second->get_avg_z());
0104 
0105 // handle rollover from pi to -pi
0106     double diffphi = phi-ntvars[0];
0107     if (diffphi > M_PI)
0108     {
0109       diffphi -= 2*M_PI;
0110     }
0111     else if (diffphi < - M_PI)
0112     {
0113       diffphi += 2*M_PI;
0114     }
0115     //  double difftheta = theta-ntvars[1];
0116     if (theta < 0)
0117       {
0118         theta += M_PI;
0119       }
0120     double difftheta = theta-ntvars[1];
0121 // theta goes from 0-PI --> no rollover problem
0122     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0123     double edep = hit_iter->second->get_edep();
0124     eall[0] += edep;
0125     detid[0] = 0;
0126     detid[1] = diffphi;
0127     detid[2] = theta-ntvars[1];
0128     detid[3] = deltasqrt;
0129     detid[4] = edep;
0130     ntup->Fill(detid);
0131     for (int i=0; i<10; i++)
0132     {
0133       if (deltasqrt < (i+1)*0.025)
0134       {
0135         emc[i]+=edep;
0136       }
0137     }
0138       }
0139     }
0140   }
0141   nodename[0] = {"G4HIT_HCALIN"};
0142   nodename[1] = {"G4HIT_ABSORBER_HCALIN"};
0143   for (int j=0; j<2;j++)
0144   {
0145     hits = findNode::getClass<PHG4HitContainer>(topNode,nodename[j]);
0146     if (hits)
0147     {
0148       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0149       for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0150       {
0151     // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0152     //           hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0153     //           hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0154     double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0155     double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0156                  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0157                 hit_iter->second->get_avg_z());
0158  // handle rollover from pi to -pi
0159     double diffphi = phi-ntvars[0];
0160     if (diffphi > M_PI)
0161     {
0162       diffphi -= 2*M_PI;
0163     }
0164     else if (diffphi < - M_PI)
0165     {
0166       diffphi += 2*M_PI;
0167     }
0168     //  double difftheta = theta-ntvars[1];
0169     if (theta < 0)
0170       {
0171         theta += M_PI;
0172       }
0173     double difftheta = theta-ntvars[1];
0174 // theta goes from 0-PI --> no rollover problem
0175     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0176     double edep = hit_iter->second->get_edep();
0177     eall[0] += edep;
0178     detid[0] = 1;
0179     detid[1] = diffphi;
0180     detid[2] = theta-ntvars[1];
0181     detid[3] = deltasqrt;
0182     detid[4] = edep;
0183     ntup->Fill(detid);
0184     for (int i=0; i<10; i++)
0185     {
0186       if (deltasqrt < (i+1)*0.025)
0187       {
0188         ih[i]+=edep;
0189       }
0190     }
0191       }
0192     }
0193   }
0194   nodename[0] = {"G4HIT_HCALOUT"};
0195   nodename[1] = {"G4HIT_ABSORBER_HCALOUT"};
0196   for (int j=0; j<2;j++)
0197   {
0198     hits = findNode::getClass<PHG4HitContainer>(topNode,nodename[j]);
0199     if (hits)
0200     {
0201       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0202       for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0203       {
0204     // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0205     //           hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0206     //           hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0207     double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0208     double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0209                  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0210                 hit_iter->second->get_avg_z());
0211 // handle rollover from pi to -pi
0212     double diffphi = phi-ntvars[0];
0213     if (diffphi > M_PI)
0214     {
0215       diffphi -= 2*M_PI;
0216     }
0217     else if (diffphi < - M_PI)
0218     {
0219       diffphi += 2*M_PI;
0220     }
0221     //  double difftheta = theta-ntvars[1];
0222     if (theta < 0)
0223       {
0224         theta += M_PI;
0225       }
0226     double difftheta = theta-ntvars[1];
0227 // theta goes from 0-PI --> no rollover problem
0228     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0229     double edep = hit_iter->second->get_edep();
0230     eall[0] += edep;
0231     detid[0] = 2;
0232     detid[1] = diffphi;
0233     detid[2] = theta-ntvars[1];
0234     detid[3] = deltasqrt;
0235     detid[4] = edep;
0236     ntup->Fill(detid);
0237     for (int i=0; i<10; i++)
0238     {
0239       if (deltasqrt < (i+1)*0.025)
0240       {
0241         oh[i]+=edep;
0242       }
0243     }
0244       }
0245     }
0246   }
0247   hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_MAGNET");
0248   if (hits)
0249   {
0250     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0251     for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0252     {
0253       // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0254       //             hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0255       //             hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0256       double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0257       double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0258                    hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0259               hit_iter->second->get_avg_z());
0260 // handle rollover from pi to -pi
0261     double diffphi = phi-ntvars[0];
0262     if (diffphi > M_PI)
0263     {
0264       diffphi -= 2*M_PI;
0265     }
0266     else if (diffphi < - M_PI)
0267     {
0268       diffphi += 2*M_PI;
0269     }
0270     //double difftheta = theta-ntvars[1];
0271     if (theta < 0)
0272       {
0273         theta += M_PI;
0274       }
0275     double difftheta = theta-ntvars[1];
0276 // theta goes from 0-PI --> no rollover problem
0277     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0278     double edep = hit_iter->second->get_edep();
0279     eall[0] += edep;
0280     detid[0] = 3;
0281     detid[1] = diffphi;
0282     detid[2] = theta-ntvars[1];
0283     detid[3] = deltasqrt;
0284     detid[4] = edep;
0285     ntup->Fill(detid);
0286       for (int i=0; i<10; i++)
0287       {
0288     if (deltasqrt < (i+1)*0.025)
0289     {
0290       mag[i]+=edep;
0291     }
0292       }
0293     }
0294   }
0295   hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_BH_1");
0296   if (hits)
0297   {
0298     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0299     for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0300     {
0301       // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0302       //             hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0303       //             hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0304       double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0305       double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0306                    hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0307               hit_iter->second->get_avg_z());
0308 // handle rollover from pi to -pi
0309     double diffphi = phi-ntvars[0];
0310     if (diffphi > M_PI)
0311     {
0312       diffphi -= 2*M_PI;
0313     }
0314     else if (diffphi < - M_PI)
0315     {
0316       diffphi += 2*M_PI;
0317     }
0318     //  double difftheta = theta-ntvars[1];
0319     if (theta < 0)
0320       {
0321         theta += M_PI;
0322       }
0323     double difftheta = theta-ntvars[1];
0324 // theta goes from 0-PI --> no rollover problem
0325     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0326     double edep = hit_iter->second->get_edep();
0327     eall[0] += edep;
0328     detid[0] = 4;
0329     detid[1] = diffphi;
0330     detid[2] = theta-ntvars[1];
0331     detid[3] = deltasqrt;
0332     detid[4] = edep;
0333     ntup->Fill(detid);
0334       for (int i=0; i<10; i++)
0335       {
0336     if (deltasqrt < (i+1)*0.025)
0337     {
0338       bh[i]+=edep;
0339     }
0340       }
0341     }
0342   }
0343   for (int i=0; i<10;i++)
0344   {
0345     float nte[6] = {0};
0346     nte[0] = i;
0347     nte[1] = emc[i];
0348     nte[2] = ih[i];
0349     nte[3] = oh[i];
0350     nte[4] = mag[i];
0351     nte[5] = bh[i];
0352     ntups->Fill(nte);
0353   }
0354   float nte[6] = {0};
0355   nte[0] = 100;
0356   nte[1] = eall[0];
0357   nte[2] = eall[1];
0358   nte[3] = eall[2];
0359   nte[4] = eall[3];
0360   nte[5] = eall[4];
0361   ntups->Fill(nte);
0362 
0363   return 0;
0364 }
0365 
0366 int
0367 ShowerSize::End(PHCompositeNode * topNode)
0368 {
0369   outfile->cd();
0370   //  ntup->Write();
0371   ntups->Write();
0372   ntupe->Write();
0373   outfile->Write();
0374   outfile->Close();
0375   delete outfile;
0376   hm->dumpHistos(_filename, "UPDATE");
0377   return 0;
0378 }
0379 
0380 void
0381 ShowerSize::AddNode(const std::string &name, const int detid)
0382 {
0383   _node_postfix.insert(name);
0384   _detid[name] = detid;
0385   return;
0386 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!