Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "HijingShowerSize.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 #include <phool/PHRandomSeed.h>
0014 
0015 #include <TFile.h>
0016 #include <TH1.h>
0017 #include <TH2.h>
0018 #include <TNtuple.h>
0019 
0020 #include <gsl/gsl_rng.h>
0021 
0022 #include<sstream>
0023 
0024 using namespace std;
0025 
0026 HijingShowerSize::HijingShowerSize(const std::string &name, const std::string &filename):
0027   SubsysReco( name ),
0028   nblocks(0),
0029   hm(nullptr),
0030   _filename(filename),
0031   ntups(nullptr),
0032   ntupe(nullptr),
0033   ntup(nullptr),
0034   outfile(nullptr)
0035 {
0036   RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
0037   seed = PHRandomSeed(); // fixed seed is handled in this funtcion
0038   cout << Name() << " random seed: " << seed << endl;
0039   gsl_rng_set(RandomGenerator,seed);
0040   return;
0041 }
0042 
0043 HijingShowerSize::~HijingShowerSize()
0044 {
0045   //  delete ntup;
0046   gsl_rng_free (RandomGenerator);
0047   delete hm;
0048 } 
0049 
0050 
0051 int
0052 HijingShowerSize::Init( PHCompositeNode* )
0053 {
0054   ostringstream hname, htit;
0055   hm = new  Fun4AllHistoManager(Name());
0056   outfile = new TFile(_filename.c_str(), "RECREATE");
0057   ntups = new TNtuple("sz"," Shower size","rad:em:hi:ho:mag:bh");
0058   ntupe = new TNtuple("truth", "The Absolute Truth", "phi:theta");
0059   ntup = new TNtuple("de", "Change in Angles", "ID:dphi:dtheta:dtotal:edep");
0060   return 0;
0061 }
0062 
0063 int
0064 HijingShowerSize::process_event( PHCompositeNode* topNode )
0065 {
0066   float ntvars[5][2] = {0};
0067   float tupvars[2];
0068   for (int i=0; i<5; i++)
0069   {
0070     ntvars[i][0] = 4*M_PI*gsl_rng_uniform_pos(RandomGenerator) - 2*M_PI;
0071     ntvars[i][1] = (1.62-1.52)*gsl_rng_uniform_pos(RandomGenerator) + 1.52;
0072     if (ntvars[i][1] < 0)
0073       {
0074     ntvars[i][1]+=M_PI;
0075       }
0076     tupvars[0] = ntvars[i][0];
0077     tupvars[1] = ntvars[i][1];
0078   ntupe->Fill(tupvars);
0079   }
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     for (int i = 0; i<5; i++)
0107     {
0108     double diffphi = phi-ntvars[i][0];
0109     if (diffphi > M_PI)
0110     {
0111       diffphi -= 2*M_PI;
0112     }
0113     else if (diffphi < - M_PI)
0114     {
0115       diffphi += 2*M_PI;
0116     }
0117     //  double difftheta = theta-ntvars[1];
0118     if (theta < 0)
0119       {
0120         theta += M_PI;
0121       }
0122     double difftheta = theta-ntvars[i][1];
0123 // theta goes from 0-PI --> no rollover problem
0124     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0125     if (deltasqrt>0.5)
0126     {
0127       continue;
0128     }
0129     double edep = hit_iter->second->get_edep();
0130     eall[0] += edep;
0131     detid[0] = 0;
0132     detid[1] = diffphi;
0133     detid[2] = difftheta;
0134     detid[3] = deltasqrt;
0135     detid[4] = edep;
0136     ntup->Fill(detid);
0137     for (int i=0; i<10; i++)
0138     {
0139       if (deltasqrt < (i+1)*0.025)
0140       {
0141         emc[i]+=edep;
0142       }
0143     }
0144     }
0145       }
0146     }
0147   }
0148   nodename[0] = {"G4HIT_HCALIN"};
0149   nodename[1] = {"G4HIT_ABSORBER_HCALIN"};
0150   for (int j=0; j<2;j++)
0151   {
0152     hits = findNode::getClass<PHG4HitContainer>(topNode,nodename[j]);
0153     if (hits)
0154     {
0155       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0156       for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0157       {
0158     // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0159     //           hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0160     //           hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0161     double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0162     double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0163                  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0164                 hit_iter->second->get_avg_z());
0165  // handle rollover from pi to -pi
0166     for (int i=0; i<5; i++)
0167     {
0168     double diffphi = phi-ntvars[i][0];
0169     if (diffphi > M_PI)
0170     {
0171       diffphi -= 2*M_PI;
0172     }
0173     else if (diffphi < - M_PI)
0174     {
0175       diffphi += 2*M_PI;
0176     }
0177     //  double difftheta = theta-ntvars[1];
0178     if (theta < 0)
0179       {
0180         theta += M_PI;
0181       }
0182     double difftheta = theta-ntvars[i][1];
0183 // theta goes from 0-PI --> no rollover problem
0184     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0185     if (deltasqrt>0.5)
0186     {
0187       continue;
0188     }
0189     double edep = hit_iter->second->get_edep();
0190     eall[0] += edep;
0191     detid[0] = 1;
0192     detid[1] = diffphi;
0193     detid[2] = difftheta;
0194     detid[3] = deltasqrt;
0195     detid[4] = edep;
0196     ntup->Fill(detid);
0197     for (int i=0; i<10; i++)
0198     {
0199       if (deltasqrt < (i+1)*0.025)
0200       {
0201         ih[i]+=edep;
0202       }
0203     }
0204     }
0205       }
0206     }
0207   }
0208   nodename[0] = {"G4HIT_HCALOUT"};
0209   nodename[1] = {"G4HIT_ABSORBER_HCALOUT"};
0210   for (int j=0; j<2;j++)
0211   {
0212     hits = findNode::getClass<PHG4HitContainer>(topNode,nodename[j]);
0213     if (hits)
0214     {
0215       PHG4HitContainer::ConstRange hit_range = hits->getHits();
0216       for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0217       {
0218     // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0219     //           hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0220     //           hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0221     double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0222     double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0223                  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0224                 hit_iter->second->get_avg_z());
0225 // handle rollover from pi to -pi
0226     for (int i=0; i<5; i++)
0227     {
0228     double diffphi = phi-ntvars[i][0];
0229     if (diffphi > M_PI)
0230     {
0231       diffphi -= 2*M_PI;
0232     }
0233     else if (diffphi < - M_PI)
0234     {
0235       diffphi += 2*M_PI;
0236     }
0237     //  double difftheta = theta-ntvars[1];
0238     if (theta < 0)
0239       {
0240         theta += M_PI;
0241       }
0242     double difftheta = theta-ntvars[i][1];
0243 // theta goes from 0-PI --> no rollover problem
0244     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0245     if (deltasqrt>0.5)
0246     {
0247       continue;
0248     }
0249     double edep = hit_iter->second->get_edep();
0250     eall[0] += edep;
0251     detid[0] = 2;
0252     detid[1] = diffphi;
0253     detid[2] = difftheta;
0254     detid[3] = deltasqrt;
0255     detid[4] = edep;
0256     ntup->Fill(detid);
0257     for (int i=0; i<10; i++)
0258     {
0259       if (deltasqrt < (i+1)*0.025)
0260       {
0261         oh[i]+=edep;
0262       }
0263     }
0264     }
0265       }
0266     }
0267   }
0268   hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_MAGNET");
0269   if (hits)
0270   {
0271     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0272     for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0273     {
0274       // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0275       //             hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0276       //             hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0277       double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0278       double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0279                    hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0280               hit_iter->second->get_avg_z());
0281 // handle rollover from pi to -pi
0282       for (int i=0; i<5; i++)
0283       {
0284     double diffphi = phi-ntvars[i][0];
0285     if (diffphi > M_PI)
0286     {
0287       diffphi -= 2*M_PI;
0288     }
0289     else if (diffphi < - M_PI)
0290     {
0291       diffphi += 2*M_PI;
0292     }
0293     //double difftheta = theta-ntvars[1];
0294     if (theta < 0)
0295       {
0296         theta += M_PI;
0297       }
0298     double difftheta = theta-ntvars[i][1];
0299 // theta goes from 0-PI --> no rollover problem
0300     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0301     if (deltasqrt>0.5)
0302     {
0303       continue;
0304     }
0305     double edep = hit_iter->second->get_edep();
0306     eall[0] += edep;
0307     detid[0] = 3;
0308     detid[1] = diffphi;
0309     detid[2] = difftheta;
0310     detid[3] = deltasqrt;
0311     detid[4] = edep;
0312     ntup->Fill(detid);
0313       for (int i=0; i<10; i++)
0314       {
0315     if (deltasqrt < (i+1)*0.025)
0316     {
0317       mag[i]+=edep;
0318     }
0319       }
0320       }
0321     }
0322   }
0323   hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_BH_1");
0324   if (hits)
0325   {
0326     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0327     for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0328     {
0329       // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0330       //             hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0331       //             hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0332       double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0333       double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0334                    hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0335               hit_iter->second->get_avg_z());
0336 // handle rollover from pi to -pi
0337       for (int i=0; i<5; i++)
0338       {
0339     double diffphi = phi-ntvars[i][0];
0340     if (diffphi > M_PI)
0341     {
0342       diffphi -= 2*M_PI;
0343     }
0344     else if (diffphi < - M_PI)
0345     {
0346       diffphi += 2*M_PI;
0347     }
0348     //  double difftheta = theta-ntvars[1];
0349     if (theta < 0)
0350       {
0351         theta += M_PI;
0352       }
0353     double difftheta = theta-ntvars[i][1];
0354 // theta goes from 0-PI --> no rollover problem
0355     double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0356     if (deltasqrt>0.5)
0357     {
0358       continue;
0359     }
0360     double edep = hit_iter->second->get_edep();
0361     eall[0] += edep;
0362     detid[0] = 4;
0363     detid[1] = diffphi;
0364     detid[2] = difftheta;
0365     detid[3] = deltasqrt;
0366     detid[4] = edep;
0367     ntup->Fill(detid);
0368       for (int i=0; i<10; i++)
0369       {
0370     if (deltasqrt < (i+1)*0.025)
0371     {
0372       bh[i]+=edep;
0373     }
0374       }
0375       }
0376     }
0377   }
0378   for (int i=0; i<10;i++)
0379   {
0380     float nte[6] = {0};
0381     nte[0] = i;
0382     nte[1] = emc[i];
0383     nte[2] = ih[i];
0384     nte[3] = oh[i];
0385     nte[4] = mag[i];
0386     nte[5] = bh[i];
0387     ntups->Fill(nte);
0388   }
0389   float nte[6] = {0};
0390   nte[0] = 100;
0391   nte[1] = eall[0];
0392   nte[2] = eall[1];
0393   nte[3] = eall[2];
0394   nte[4] = eall[3];
0395   nte[5] = eall[4];
0396   ntups->Fill(nte);
0397 
0398   return 0;
0399 }
0400 
0401 int
0402 HijingShowerSize::End(PHCompositeNode * topNode)
0403 {
0404   outfile->cd();
0405   //  ntup->Write();
0406   ntups->Write();
0407   ntupe->Write();
0408   outfile->Write();
0409   outfile->Close();
0410   delete outfile;
0411   hm->dumpHistos(_filename, "UPDATE");
0412   return 0;
0413 }
0414 
0415 void
0416 HijingShowerSize::AddNode(const std::string &name, const int detid)
0417 {
0418   _node_postfix.insert(name);
0419   _detid[name] = detid;
0420   return;
0421 }
[ 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!