EnergyCorrelatorsJets/CalorimeterEnergyCorrelatorCommonTools/BuildMetaTowers.h

0001 #ifndef __BUILDMETATOWERS_H__
0002 #define __BUILDMETATOWERS_H__
0003
0004 //////////////////////////////////////////////////////////////////////////
0005 //                                  //
0006 //          Build Meta Towers               //
0007 //                                  //
0008 //  Does what is says on the tin, take in calorimeters, applies a   //
0009 //  vertex based correction and combines the towers to make a   //
0010 //  single "shifted em+hadronic calorimeter" w/ HCAL bins       //
0011 //  Works with vandy skimmer data and Tower info            //
0012 //                                  //
0013 //  Authors:    Ben Kimelman, Skaydi                //
0014 //  First commit:   9 March 2026                    //
0015 //  This commit:    10 March 2026                   //
0016 //  version:    v1.0                        //
0017 //                                  //
0018 //  Notes on this commit: First commit              //
0019 //                                  //
0020 //////////////////////////////////////////////////////////////////////////
0021
0022
0023 //c++ classes
0024
0025 #include <math.h>
0026 #include <vector>
0027 #include <array>
0028 #include <string>
0029
0030 //Vandy classes
0031 #include <vandyclasses/Tower.h>
0032 #include <vadyclasses/EventInfo.h>
0033
0034 //Calo classes
0035 #include <calobase/TowerInfoContainer.h>
0036 #include <calobase/TowerInfoContainerv1.h>
0037 #include <calobase/TowerInfoContainerv2.h>
0038 #include <calobase/TowerInfov2.h>
0039 #include <calobase/TowerInfov1.h>
0040 #include <calobase/TowerInfo.h>
0041 #include <calobase/RawTowerDefs.h>
0042 #include <calobase/RawTowerContainer.h>
0043 #include <calobase/RawTowerGeomContainer.h>
0044 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0045
0046 //G4 objects
0047 #include <g4main/PHG4Particle.h>
0048 #include <g4main/PHG4Hit.h>
0049 #include <g4main/PHG4TruthInfoContainer.h>
0050
0051 //phool
0052 #include <phool/PHCompositeNode.h>
0053 #include <phool/getClass.h>
0054 #include <phool/PHNodeOperation.h>
0055
0056 class BuildMetaTowers
0057 {
0058     public:
0059         BuildMetaTowers(
0060                 float radius = 1.245 /*average inner hcal radius*/,
0061                 std::string input="VandyClass" /*Vandy Class, TowerInfo or Array*/
0062                    )
0063         {
0064             //assume that if radius given is small it is given in m, else
0065             if(radius > 50 ) this->R = radius;
0066             if(radius < 50 ) this->R = radius * 100;
0067             this->T = input;
0068         }
0069         ~BuildMetaTowers(){};
0070
0071         //struct for Array input
0072         struct TowerArrayEntry
0073         {
0074             double Energy;
0075             double phi; //center of tower phi
0076             double eta; //center of tower eta
0077         };
0078         enum CALO
0079         {
0080             META    = 0,
0081             EMCAL   = 1,
0082             IHCAL   = 2,
0083             OHCAL   = 3
0084         };
0085         //Common Methods
0086         void RunMetaTowerBuilder(float zVTX)
0087         {
0088             for( int i = 0; i < 24; i++ ) shiftedetaEdges.at(i)=calculateEtaShift(etaEdges.at(i), zVTX);
0089             for( int i = 0; i < 24; i++ )
0090             {
0091                 double eta =10.;
0092                 if(i < (int) shiftedetaEdges.size() - 1 )
0093                     eta = 0.5*(shiftedEtaEdges.at(i),  shiftedEtaEdges.at(i+1));
0094                 else
0095                     eta = 0.5*(shiftedEtaEdges.at(i),  1.1);
0096                 for( int j = 0; j < 64; j++)
0097                 {
0098                     double pt =100.;
0099                     TowerArrayEntry* tower =
0100                         new TowerArrayEntry* {0., shiftedEtaEdges.at(i), phiEdges.at(i)};
0101                     if(i < (int) phiEdges.size() - 1 )
0102                         phi = 0.5*(phiEdges.at(i),  phiEdges.at(i+1));
0103                     else
0104                         phi = 0.5*(phiEdges.at(i),  phiEdges.at(0));
0105                     if(phi > M_PI)  phi += -2*M_PI;
0106                     if(phi < -M_PI) phi +=  2*M_PI;
0107                     tower->eta = eta;
0108                     tower->phi = phi;
0109                     int index  = CalculateIndex(tower);
0110                     MetaTower->at(index)=tower;
0111                 }
0112
0113             }
0114             for(int i=0; i<(int)EMReTowers->size(); i++) AddMetaTower(EMReTowers->at(i));
0115             for(int i=0; i<(int)IHCaTowers->size(); i++) AddMetaTower(IHCaTowers->at(i));
0116             for(int i=0; i<(int)OHCaTowers->size(); i++) AddMetaTower(OHCaTowers->at(i));
0117             return;
0118
0119         }
0120         //Fun4All load in
0121         void GetFun4AllTowers(
0122                 TowerInfoContainerv2 CAL,
0123                 RawTowerGeomContainer_Cylinderv1 geom,
0124                 std::array<TowerArrayEntry*, 1536>* Towers
0125                 double zVTX,
0126                 CALO calo
0127                 )
0128         {
0129             for(int t = 0; t < (int)CAL.size(); t++)
0130             {
0131                 auto key = CAL.encode_key(n) ;
0132                 auto f4tower    = CAL.get_tower_at_channel(t);
0133                 int phibin  = CAL.getTowerPhiBin(key);
0134                 int etabin  = CAL.getTowerEtabin(key);
0135                 double phi  = geom.get_phicenter(phibin);
0136                 double eta  = geom.get_etaecenter(etabin);
0137                 double E    = f4tower->get_energy();
0138                 eta         = CalculateEtaShift(eta, zVTX, calo);
0139                 TowerArrayEntry* tower =
0140                     new TowerArrayEntry {E, phi, eta};
0141                 int index   = CalculateIndex(tower, true);
0142                 Towers->at(index) = tower;
0143             }
0144             return;
0145         }
0146         void GetEMCALTowers(
0147                 TowerInfoContainerv2 EMCAL,
0148                 RawTowerGeomContainer_Cylinderv1 EMCAL_geom,
0149                 bool isRetower=false,
0150                 double zVTX = 0.
0151                 ) //if using Fun4AllTower objects
0152         {
0153             if(!isRetower){
0154                 std::cout<<"Please use the retowered EMCAL" <<std::endl;
0155                 return;
0156             }
0157             CALO cal = CALO::EMCAL;
0158             GetFun4AllTowers(EMCAL, EMCAL_geom, EMReTowers, zVTX, cal);
0159             return;
0160         }
0161         void GetHCALTowers(
0162                 TowerInfoContainerv2 HCAL,
0163                 RawTowerGeomContainer_Cylinderv1 HCAL_geom,
0164                 bool outer=false,
0165                 double zVTX = 0.
0166                 )
0167         {
0168             std::array<TowerEntryArray*, 1536>* HCaTowers = NULL;
0169             CALO cal = CALO::OHCAL;
0170             if(outer){
0171                 HCaTowers   = OHCaTowers;
0172                 cal         = CALO::OHCAL;
0173             }
0174             else{
0175                 HCaTowers   = IHCaTowers;
0176                 cal         = CALO::IHCAL;
0177             }
0178             GetFun4AllTowers(HCAL, HCAL_geom, HCaTowers, zVTX, cal);
0179             return;
0180         }
0181
0182             //if using Fun4AllTower objects
0183         //Vandy Class load in
0184         void GetEMCALTowers(
0185                 std::vector<Tower*> EMCAL,
0186                 bool isRetower=false,
0187                 double zVTX = 0.
0188                 )
0189         {
0190             if(!isRetower){
0191                 std::cout<<"Please use retowered"<<std::endl;
0192                 return;
0193             }
0194             CALO calo   = CALO::EMCAL;
0195             GetVandyTowers(EMCAL, EMReTowers, zVTX, calo);
0196             return;
0197         }
0198         void GetHCALTowers(
0199                 std::vector<Tower*> HCAL,
0200                 bool outer=false,
0201                 double zVTX = 0.
0202                 );
0203         {
0204             if(outer)
0205             {
0206                 CALO calo   = CALO::OHCAL;
0207                 GetVandyTowers(HCAL, OHCaTowers, zVTX, calo);
0208             }
0209             else if(!outer)
0210             {
0211                 CALO calo   = CALO::IHCAL;
0212                 GetVandyTowers(HCAL, IHCaTowers, zVTX, calo);
0213             }
0214             return;
0215         }
0216         void GetVandyTowers(
0217                 std::vector<Tower*> CAL,
0218                 std::array<TowerArrayEntry*, 1536>* output_array
0219                 double zVTX,
0220                 CALO calo
0221                 )
0222         {
0223             for(int i = 0; i<(int)CAL.size(); i++)
0224             {
0225                 Tower* tw   = CAL.at(i);
0226                 double E    = tw->e();
0227                 double px   = tw->px();
0228                 double py   = tw->py();
0229                 double pz   = tw->pz();
0230                 double p    = std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2));
0231                 double phi  = std::atan2(px, py);
0232                 double eta  = std::atanh(pz/p);
0233                 eta         = CalculateEtaShift(eta, zVTX, calo);
0234                 TowerArrayEntry* tower = new TowerArrayEntry {E, phi, eta};
0235                 int index   = CalculateIndex(tower, true);
0236                 output_array->at(i) = tower;
0237             }
0238             return;
0239         }
0240
0241         //Just basic array
0242         void GetEMCALTowers(
0243                 std::vector<TowerArrayEntry*> EMCAL,
0244                     bool isRetower = false,
0245                 double zVTX = 0.
0246                 )
0247         {
0248             CALO calo = CALO::EMCAL;
0249             for(int i = 0; i <(int)EMCAL->size(); i++)
0250             {
0251                 EMReTowers->at(i)   = EMCAL->at(i);
0252                 EMReTowers->at(i)->eta  = calculateEtaShift(EMReTower->at(i)->eta, zVTX, calo);
0253             }
0254             return;
0255         }
0256
0257         void GetHCALTowers(
0258                 std::vector<TowerArrayEntry*> HCAL,
0259                 bool outer=false,
0260                 double zVTX = 0.
0261                 )
0262         {
0263             CALO calo = CALO::OHCAL;
0264             std::array<TowerArrayEntry*>* HC = OHCaTowers;
0265             if(!outer){
0266                     calo    = CALO::IHCAL;
0267                 HC  = IHCaTowers;
0268             for(int i = 0; i <(int)HCAL->size(); i++)
0269             {
0270                 HC->at(i)->eta  = calculateEtaShift(HCAL->at(i)->eta, zVTX, calo);
0271             }
0272             return;
0273         }
0274         //just returns configs
0275         void setRadius  ( float r = 1.245 ) { this->R = r ; } ;
0276             float getRadius ( ) { return this->R; };
0277
0278         void        setInput
0279                 ( std::string input = "VandyClass" ) { this->T = input; };
0280         std::string     getInput() { return this->T; };
0281
0282         //return the final object
0283         std::array<TowerArrayEntry*, 1536>* getMetaTowers() { return this->MetaTowers; };
0284
0285     private:
0286         //Private variables
0287         float       R { 1245 }; //IHCAL half radius in cm
0288         std::string     T { "VandyClass" }; //What input Type to expect
0289         const float     R_OHCAL_Outer {269};
0290         const float     R_OHCAL_Inner {182};
0291         const float     R_IHCAL_Outer {137};
0292         const float     R_IHCAL_Inner {116};
0293         const float     R_EMCAL_Outer {116};
0294         const float R_EMCAL_Inner { 90};
0295         const float R_OHCAL_MID { 1/2.*(R_OHCAL_Outer + R_OHCAL_Inner) };
0296         const float R_IHCAL_MID { 1/2.*(R_IHCAL_Outer + R_IHCAL_Inner) };
0297         const float R_EMCAL_MID { 1/2.*(R_EMCAL_Outer + R_EMCAL_Inner) };
0298         std::array <TowerArrayEntry*, 1536>* MetaTowers = new std::array <TowerArrayEntry*, 1536> {};
0299         std::array <TowerArrayEntry*, 1536>* EMReTowers = new std::array <TowerArrayEntry*, 1536> {};
0300         std::array <TowerArrayEntry*, 1536>* IHCaTowers = new std::array <TowerArrayEntry*, 1536> {};
0301         std::array <TowerArrayEntry*, 1536>* OHCaTowers = new std::array <TowerArrayEntry*, 1536> {};
0302
0303         //HCAL eta-phi physical geometry
0304         const std::array <double, 24> etaEdges
0305             { -1.1000000, -1.0083333, -0.91666667, -0.82500000, -0.73333333, -0.64166667,
0306                 -0.55000000, -0.45833333, -0.36666667, -0.27500000, -0.18333333, -0.091666667,
0307                 1.1102230e-16, 0.091666667, 0.18333333, 0.27500000, 0.36666667, 0.45833333,
0308                 0.55000000, 0.64166667, 0.73333333, 0.82500000, 0.91666667, 1.0083333, 1.1000000 };
0309
0310         const std::array <double, 64> phiEdges
0311             { -0.053619781, 0.044554989, 0.14272976, 0.24090453, 0.33907930, 0.43725407,
0312                 0.53542884, 0.63360361, 0.73177838, 0.82995315, 0.92812792, 1.0263027,
0313                 1.1244775, 1.2226522, 1.3208270, 1.4190018, 1.5171765, 1.6153513, 1.7135261,
0314                 1.8117009, 1.9098756, 2.0080504, 2.1062252, 2.2043999, 2.3025747, 2.4007495,
0315                 2.4989242, 2.5970990, 2.6952738, 2.7934486, 2.8916233, 2.9897981, 3.0879729,
0316                 3.1861476, 3.2843224, 3.3824972, 3.4806720, 3.5788467, 3.6770215, 3.7751963,
0317                 3.8733710, 3.9715458, 4.0697206, 4.1678953, 4.2660701, 4.3642449, 4.4624197,
0318                 4.5605944, 4.6587692, 4.7569440, 4.8551187, 4.9532935, 5.0514683, 5.1496431,
0319                 5.2478178, 5.3459926, 5.4441674, 5.5423421, 5.6405169, 5.7386917, 5.8368664,
0320                 5.9350412, 6.0332160, 6.1313908, 6.2295655 };
0321
0322         std::array<double ,24> shiftedetaEdges {};
0323
0324         //Private Methods
0325         double calculateEtaShift(double eta, double zVtx)
0326         {
0327             double z    = R*std::sinh(eta);
0328             double zshift   = z-zVtx;
0329             zshift      = zshift / R;
0330             double etashift = std::asinh(zshift);
0331             return etashift;
0332         }
0333         double calculateEtaShift(double eta, double zVtx, CAlO calo)
0334         {
0335             double r = this->R;
0336             if(         calo == CALO::EMCAL )   r = R_EMCAL_MID;
0337             else if(    calo == CALO::IHCAL )   r = R_IHCAL_MID;
0338             else if(    calo == CALO::OHCAL )   r = R_OHCAL_MID;
0339             else                    r = this->R;
0340             double z    = R*std::sinh(eta);
0341             double zshift   = z-zVtx;
0342             zshift      = zshift / R;
0343             double etashift = std::asinh(zshift);
0344             return etashift;
0345         }
0346
0347
0348         int calculateIndex(TowerArrayEntry tower, bool useSlant)
0349         {
0350             int etabin  = 0;
0351             int phibin  = 0;
0352             etabin      = LookupBin(tower.eta, useSlant, true);
0353             phibin      = LookupBin(tower.phi, useSlant, false);
0354             int index = etabin * 64 + phibin ;
0355             return index;
0356         }
0357         int LookupBin(double towervalue, bool useSlant, bool isEta)
0358         {
0359             std::vector<double> tempBottom;
0360             int binN=0;
0361             if(isEta && useSlant)
0362             {
0363                 for(int i = 0; i<(int) shiftedEtaEdges.size(); i++){
0364                     tempBottom.push_back(shiftedEtaEdges.at(i));
0365                 }
0366             }
0367             else if (isEta)
0368             {
0369                 for(int i = 0; i<(int) EtaEdges.size(); i++){
0370                     tempBottom.push_back(EtaEdges.at(i));
0371                 }
0372             }
0373             else
0374             {
0375                 for(int i = 0; i<(int) PhiEdges.size(); i++){
0376                     tempBottom.push_back(PhiEdges.at(i));
0377                 }
0378             }
0379             for(int i = 0; i<(int)tempBottom.size(); i++)
0380             {
0381                 if(towervalue >= tempBottom.at(i)){
0382                     binN = i;
0383                     continue;
0384                 }
0385                 if(towervalue < tempBottom.at(i)){
0386                     break;
0387                 }
0388             }
0389             return binN;
0390         }
0391         void decodeIndex(int index, int* ebr, int* pbr)
0392         {
0393
0394             int etabin  = index / 64;
0395             int phibin  = index % 64;
0396             ebr = &etabin;
0397             pbr = &phibin;
0398             return;
0399         }
0400
0401         void shiftBinMins(double zVtx)
0402         {
0403             for(int i= 0; i<(int) etaEdges.size(); i++)
0404             {
0405                 shiftedetaEdges[i]=calculateEtaShift(etaEdges[i], zVtx);
0406             }
0407             return;
0408         }
0409         void shiftEMCAL(std::vector<TowerArrayEntry*> emcal, bool isRetower, float zVTX)
0410         {
0411             if(!isRetower || (int) emcal.size() > 1536)
0412             {
0413                 std::vector<TowerArrayEntry*> rtEM = RetowerEMCAL(emcal);
0414                 emcal.clear();
0415                 emcal = rtEM;
0416             }
0417             for(int i = 0; i<(int) emcal.size(); i++)
0418             {
0419                 TowerArrayEntry* nTower = shiftTower(emcal.at(i), CALO::EMC, zVTX);
0420                 int N = calculateIndex( *nTower );
0421                 EMReTowers->at(N) = nTower;
0422             }
0423             return;
0424         }
0425
0426         void shiftIHCAL(std::vector<TowerArrayEntry*> ihcal, float zVTX);
0427         {
0428             for(int i = 0; i<(int) ihcal.size(); i++)
0429             {
0430                 TowerArrayEntry* nTower = shiftTower(ihcal.at(i), CALO::IHC, zVTX);
0431                 int N = calculateIndex( *nTower );
0432                 IHCaTowers->at(N) = nTower;
0433             }
0434             return;
0435         }
0436         void shiftOHCAL(std::vector<TowerArrayEntry*> ohcal, float zVTX);
0437         {
0438             for(int i = 0; i<(int) ohcal.size(); i++)
0439             {
0440                 TowerArrayEntry* nTower = shiftTower(ohcal.at(i), CALO::OHC, zVTX);
0441                 int N = calculateIndex( *nTower );
0442                 OHCaTowers->at(N) = nTower;
0443             }
0444             return;
0445         }
0446         void addMetaTower(TowerArrayEntry tower)
0447         {
0448             double eta_val  = tower.eta;
0449             double phi_val  = tower.phi;
0450             int etaBin  = findEtaBin(eta_val);
0451             int phiBin  = findPhiBin(phi_val);
0452             int N       = calculateIndex(tower);
0453             if (!MetaTowers->at(N)) MetaTowers->at(N) =
0454                 new TowerArrayEntry* { 0, getPhiCenter(phiBin), getEtaCenter(etaBin) };
0455             MetaTowers->at(N)->E += tower.E;
0456             return;
0457         }
0458         TowerArrayEntry* shiftTower(TowerArrayEntry* tower, CALO calo, float zVTX)
0459         {
0460
0461             tower->eta = CalculateEtaShift(tower->eta, zVTX, calo);
0462             return tower;
0463         }
0464
0465 };
0466
0467 #endif