Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-09-28 07:03:00 

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