Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-18 09:11:04 

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