Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-07-18 08:51:37 

0001 // Copyright 2023, Friederike Bock
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003 //
0004 //  Sections Copyright (C) 2023 Friederike Bock
0005 //  under SPDX-License-Identifier: LGPL-3.0-or-later
0006 
0007 #include <vector>
0008 #include <TVector3.h>
0009 
0010 struct towersStrct {
0011   towersStrct()
0012       : energy(0)
0013       , time(0)
0014       , posx(0)
0015       , posy(0)
0016       , posz(0)
0017       , cellID(0)
0018       , cellIDx(-1)
0019       , cellIDy(-1)
0020       , cellIDz(-1)
0021       , tower_trueID(-10000)
0022       , tower_clusterIDA(-1)
0023       , tower_clusterIDB(-1) {}
0024   float energy;
0025   float time;
0026   float posx;
0027   float posy;
0028   float posz;
0029   int cellID;
0030   int cellIDx;
0031   int cellIDy;
0032   int cellIDz;
0033   int tower_trueID;
0034   int tower_clusterIDA;
0035   int tower_clusterIDB;
0036 };
0037 
0038 bool acompare(towersStrct lhs, towersStrct rhs) { return lhs.energy > rhs.energy; }
0039 
0040 struct clustersStrct {
0041   clustersStrct()
0042       : cluster_E(0.)
0043       , cluster_seed(0.)
0044       , cluster_Eta(-10.)
0045       , cluster_Phi(-10.)
0046       , cluster_X(0.)
0047       , cluster_Y(0.)
0048       , cluster_Z(0.)
0049       , cluster_M02(0.)
0050       , cluster_M20(0.)
0051       , cluster_NTowers(0)
0052       , cluster_trueID(-10000)
0053       , cluster_NtrueID(0) {}
0054   float cluster_E;
0055   float cluster_seed;
0056   float cluster_Eta;
0057   float cluster_Phi;
0058   float cluster_X;
0059   float cluster_Y;
0060   float cluster_Z;
0061   float cluster_M02;
0062   float cluster_M20;
0063   int cluster_NTowers;
0064   int cluster_trueID;
0065   int cluster_NtrueID;
0066   std::vector<towersStrct> cluster_towers;
0067 };
0068 
0069 bool acompareCl(clustersStrct lhs, clustersStrct rhs) { return lhs.cluster_E > rhs.cluster_E; }
0070 
0071 //**************************************************************************************************************//
0072 //**************************************************************************************************************//
0073 // find clusters with common edges or corners, separate if energy increases in neighboring cell
0074 //**************************************************************************************************************//
0075 //**************************************************************************************************************//
0076 clustersStrct findMACluster(
0077     float seed,                                    // minimum seed energy
0078     float /* agg */,                               // minimum aggregation energy
0079     std::vector<towersStrct>& input_towers_temp,   // temporary full tower array
0080     std::vector<towersStrct>& cluster_towers_temp, // towers associated to cluster
0081     //                               std::vector<int> clslabels_temp              // MC labels in cluster
0082     float aggMargin = 1.0 // aggregation margin
0083 ) {
0084   clustersStrct tempstructC;
0085   if (input_towers_temp.at(0).energy > seed) {
0086     //     std::cout << "new cluster" << std::endl;
0087     // fill seed cell information into current cluster
0088     tempstructC.cluster_E       = input_towers_temp.at(0).energy;
0089     tempstructC.cluster_seed    = input_towers_temp.at(0).energy;
0090     tempstructC.cluster_NTowers = 1;
0091     tempstructC.cluster_NtrueID = 1;
0092     tempstructC.cluster_trueID  = input_towers_temp.at(0).tower_trueID; // TODO save all MC labels?
0093     cluster_towers_temp.push_back(input_towers_temp.at(0));
0094     //     clslabels_temp.push_back(input_towers_temp.at(0).tower_trueID);
0095     //     std::cout  << "seed: "<<  input_towers_temp.at(0).cellIDx << "\t" << input_towers_temp.at(0).cellIDy
0096     //                   << "\t" << input_towers_temp.at(0).cellIDz << "\t E:"<< tempstructC.cluster_E << std::endl;
0097 
0098     // remove seed tower from sample
0099     input_towers_temp.erase(input_towers_temp.begin());
0100     for (int tit = 0; tit < (int)cluster_towers_temp.size(); tit++) {
0101       // Now go recursively to all neighbours and add them to the cluster if they fulfill the conditions
0102       int iEtaTwr = cluster_towers_temp.at(tit).cellIDx;
0103       int iPhiTwr = cluster_towers_temp.at(tit).cellIDy;
0104       int iLTwr   = cluster_towers_temp.at(tit).cellIDz;
0105       for (int ait = 0; ait < (int)input_towers_temp.size(); ait++) {
0106         int iEtaTwrAgg = input_towers_temp.at(ait).cellIDx;
0107         int iPhiTwrAgg = input_towers_temp.at(ait).cellIDy;
0108         int iLTwrAgg   = input_towers_temp.at(ait).cellIDz;
0109 
0110         int deltaL    = TMath::Abs(iLTwrAgg - iLTwr);
0111         int deltaPhi  = TMath::Abs(iPhiTwrAgg - iPhiTwr);
0112         int deltaEta  = TMath::Abs(iEtaTwrAgg - iEtaTwr);
0113         bool neighbor = (deltaL + deltaPhi + deltaEta == 1);
0114         bool corner2D = (deltaL == 0 && deltaPhi == 1 && deltaEta == 1) ||
0115                         (deltaL == 1 && deltaPhi == 0 && deltaEta == 1) ||
0116                         (deltaL == 1 && deltaPhi == 1 && deltaEta == 0);
0117         //         first condition asks for V3-like neighbors, while second condition also checks diagonally attached towers
0118         if (neighbor || corner2D) {
0119 
0120           // only aggregate towers with lower energy than current tower
0121 
0122           if (input_towers_temp.at(ait).energy >= (cluster_towers_temp.at(tit).energy + aggMargin))
0123             continue;
0124           tempstructC.cluster_E += input_towers_temp.at(ait).energy;
0125           tempstructC.cluster_NTowers++;
0126           cluster_towers_temp.push_back(input_towers_temp.at(ait));
0127           //           if(!(std::find(clslabels_temp.begin(), clslabels_temp.end(), input_towers_temp.at(ait).tower_trueID) != clslabels_temp.end())){
0128           //             tempstructC.cluster_NtrueID++;
0129           //             clslabels_temp.push_back(input_towers_temp.at(ait).tower_trueID);
0130           //           }
0131           //           std::cout << "aggregated: "<< iEtaTwrAgg << "\t" << iPhiTwrAgg << "\t" << iLTwrAgg << "\t E:" << input_towers_temp.at(ait).energy << "\t reference: "<< refC << "\t"<< iEtaTwr << "\t" << iPhiTwr << "\t" << iLTwr << "\t cond.: \t"<< neighbor << "\t" << corner2D << "\t  diffs: " << deltaEta << "\t" << deltaPhi << "\t" << deltaL<< std::endl;
0132 
0133           input_towers_temp.erase(input_towers_temp.begin() + ait);
0134           ait--;
0135         }
0136       }
0137     }
0138   }
0139   return tempstructC;
0140 }
0141 
0142 // ANCHOR function to determine shower shape
0143 float* CalculateM02andWeightedPosition(std::vector<towersStrct> cluster_towers,
0144                                        float cluster_E_calc, float weight0) {
0145   static float returnVariables[8]; //0:M02, 1:M20, 2:eta, 3: phi
0146   float w_tot = 0;
0147   std::vector<float> w_i;
0148   TVector3 vecTwr;
0149   TVector3 vecTwrTmp;
0150   float w_0 = weight0;
0151 
0152   vecTwr = {0., 0., 0.};
0153   //calculation of weights and weighted position vector
0154   for (int cellI = 0; cellI < (int)cluster_towers.size(); cellI++) {
0155     w_i.push_back(TMath::Max(
0156         (float)0, (float)(w_0 + TMath::Log(cluster_towers.at(cellI).energy / cluster_E_calc))));
0157     w_tot += w_i.at(cellI);
0158     if (w_i.at(cellI) > 0) {
0159       vecTwrTmp = TVector3(cluster_towers.at(cellI).posx, cluster_towers.at(cellI).posy,
0160                            cluster_towers.at(cellI).posz);
0161       vecTwr += w_i.at(cellI) * vecTwrTmp;
0162     }
0163   }
0164   // correct Eta position for average shift in calo
0165   returnVariables[2] = vecTwr.Eta();
0166   returnVariables[3] =
0167       vecTwr.Phi(); //(vecTwr.Phi()<0 ? vecTwr.Phi()+TMath::Pi() : vecTwr.Phi()-TMath::Pi());
0168   vecTwr *= 1. / w_tot;
0169   //     std::cout << "Cluster: X: "<< vecTwr.X() << "\t" << " Y: "<< vecTwr.Y() << "\t" << " Z: "<< vecTwr.Z() << std::endl;
0170   returnVariables[4] = vecTwr.X();
0171   returnVariables[5] = vecTwr.Y();
0172   returnVariables[6] = vecTwr.Z();
0173 
0174   //calculation of M02
0175   float delta_phi_phi[4] = {0};
0176   float delta_eta_eta[4] = {0};
0177   float delta_eta_phi[4] = {0};
0178   float dispersion       = 0;
0179 
0180   for (int cellI = 0; cellI < (int)cluster_towers.size(); cellI++) {
0181     int iphi = cluster_towers.at(cellI).cellIDy;
0182     int ieta = cluster_towers.at(cellI).cellIDx;
0183     delta_phi_phi[1] += (w_i.at(cellI) * iphi * iphi) / w_tot;
0184     delta_phi_phi[2] += (w_i.at(cellI) * iphi) / w_tot;
0185     delta_phi_phi[3] += (w_i.at(cellI) * iphi) / w_tot;
0186 
0187     delta_eta_eta[1] += (w_i.at(cellI) * ieta * ieta) / w_tot;
0188     delta_eta_eta[2] += (w_i.at(cellI) * ieta) / w_tot;
0189     delta_eta_eta[3] += (w_i.at(cellI) * ieta) / w_tot;
0190 
0191     delta_eta_phi[1] += (w_i.at(cellI) * ieta * iphi) / w_tot;
0192     delta_eta_phi[2] += (w_i.at(cellI) * iphi) / w_tot;
0193     delta_eta_phi[3] += (w_i.at(cellI) * ieta) / w_tot;
0194 
0195     vecTwrTmp = TVector3(cluster_towers.at(cellI).posx, cluster_towers.at(cellI).posy,
0196                          cluster_towers.at(cellI).posz);
0197     // scale cluster position to z-plane
0198     vecTwr *= std::abs(vecTwrTmp.Z() / vecTwr.Z());
0199     float dx2 = pow(vecTwrTmp.X() - vecTwr.X(), 2);
0200     float dy2 = pow(vecTwrTmp.Y() - vecTwr.Y(), 2);
0201     float dz2 = pow(vecTwrTmp.Z() - vecTwr.Z(), 2);
0202     dispersion += (w_i.at(cellI) * (dx2 + dy2 + dz2)) / w_tot;
0203   }
0204   returnVariables[7] = dispersion;
0205   delta_phi_phi[0]   = delta_phi_phi[1] - (delta_phi_phi[2] * delta_phi_phi[3]);
0206   delta_eta_eta[0]   = delta_eta_eta[1] - (delta_eta_eta[2] * delta_eta_eta[3]);
0207   delta_eta_phi[0]   = delta_eta_phi[1] - (delta_eta_phi[2] * delta_eta_phi[3]);
0208 
0209   float calcM02 = 0.5 * (delta_phi_phi[0] + delta_eta_eta[0]) +
0210                   TMath::Sqrt(0.25 * TMath::Power((delta_phi_phi[0] - delta_eta_eta[0]), 2) +
0211                               TMath::Power(delta_eta_phi[0], 2));
0212   float calcM20 = 0.5 * (delta_phi_phi[0] + delta_eta_eta[0]) -
0213                   TMath::Sqrt(0.25 * TMath::Power((delta_phi_phi[0] - delta_eta_eta[0]), 2) +
0214                               TMath::Power(delta_eta_phi[0], 2));
0215   //     std::cout << "M02_calc: " << calcM02 << "\t\t = 0.5 * ( " << delta_phi_phi[0] <<" + "<<delta_eta_eta[0]<<" ) + TMath::Sqrt( 0.25 * TMath::Power( ( "<<delta_phi_phi[0]<<" - "<<delta_eta_eta[0]<<" ), 2 ) + TMath::Power( "<<delta_eta_phi[0]<<", 2 ) ) "<< std::endl;
0216   returnVariables[0] = calcM02;
0217   returnVariables[1] = calcM20;
0218   return returnVariables;
0219 }
[ 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!