ATLAS Offline Software
Loading...
Searching...
No Matches
CaloGeometryLookup.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3*/
4
7#include <TTree.h>
8#include <TVector2.h>
9#include <TCanvas.h>
10#include <TH2D.h>
11#include <TGraphErrors.h>
12#include <TVector3.h>
13#include <TLegend.h>
14
15#include "CaloDetDescr/CaloDetDescrElement.h"
16//#include "ISF_FastCaloSimParametrization/CaloDetDescrElement.h"
17#include "CaloGeoHelpers/CaloSampling.h"
19//#include "TMVA/Tools.h"
20//#include "TMVA/Factory.h"
21
22using namespace std;
23
25{
26 m_mineta=+10000;
27 m_maxeta=-10000;
28 m_minphi=+10000;
29 m_maxphi=-10000;
30
31 m_mineta_raw=+10000;
32 m_maxeta_raw=-10000;
33 m_minphi_raw=+10000;
34 m_maxphi_raw=-10000;
35
36 m_mindeta=10000;
37 m_mindphi=10000;
38
43
46 m_deta_double =0.;
47 m_dphi_double =0.;
48}
49
51= default;
52
54{
55 if(m_cells.empty()) return false;
56 for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
57 const CaloDetDescrElement* cell=ic->second;
58 if(ref->IsCompatible(cell)) return true;
59 }
60 return false;
61}
62
64{
65 for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
66 const CaloDetDescrElement* cell=ic->second;
67 ref->add(cell);
68 }
69}
70
71
73{
74 if(m_cells.empty()) return true;
75 t_cellmap::iterator ic=m_cells.begin();
76 const CaloDetDescrElement* refcell=ic->second;
77 int sampling=refcell->getSampling();
78 if(cell->getSampling()!=sampling) return false;
79 if(cell->eta_raw()*refcell->eta_raw()<0) return false;
80 if(sampling<21) { // keep only cells reasonable close in eta to each other
81 if(TMath::Abs(cell->deta()-refcell->deta())>0.001) return false;
82 if(TMath::Abs(cell->dphi()-refcell->dphi())>0.001) return false;
83
84 if( cell->eta()<mineta()-2.1*cell->deta() ) return false;
85 if( cell->eta()>maxeta()+2.1*cell->deta() ) return false;
86
87 double delta_eta=(cell->eta_raw()-refcell->eta_raw())/cell->deta();
88 double delta_phi=(cell->phi_raw()-refcell->phi_raw())/cell->dphi();
89 if(TMath::Abs(delta_eta-TMath::Nint(delta_eta)) > 0.01) return false;
90 if(TMath::Abs(delta_phi-TMath::Nint(delta_phi)) > 0.01) return false;
91
92 /*
93 cout<<"Compatible: s="<<cell->getSampling()<<"~"<<refcell->getSampling()<<"; ";
94 cout<<"eta="<<cell->eta_raw()<<","<<refcell->eta_raw()<<"; ";
95 cout<<"phi="<<cell->phi_raw()<<","<<refcell->phi_raw()<<"; ";
96 cout<<"deta="<<cell->deta()<<"~"<<refcell->deta()<<"; ";
97 cout<<"dphi="<<cell->dphi()<<"~"<<refcell->dphi()<<";";
98 cout<<"mineta="<<mineta()<<", maxeta="<<maxeta()<<"; ";
99 cout<<endl;
100 */
101 } else {
102 //FCal is not sufficiently regular to use submaps with regular mapping
103 return true;
104 /*
105 if(TMath::Abs(cell->dx()-refcell->dx())>0.001) return false;
106 if(TMath::Abs(cell->dy()-refcell->dy())>0.001) return false;
107 if( cell->x()<minx()-20*cell->dx() ) return false;
108 if( cell->x()>maxx()+20*cell->dx() ) return false;
109 if( cell->y()<miny()-20*cell->dy() ) return false;
110 if( cell->y()>maxy()+20*cell->dy() ) return false;
111
112 double delta_x=(cell->x_raw()-refcell->x_raw())/cell->dx();
113 double delta_y=(cell->y_raw()-refcell->y_raw())/cell->dy();
114 if(TMath::Abs(delta_x-TMath::Nint(delta_x)) > 0.01) return false;
115 if(TMath::Abs(delta_y-TMath::Nint(delta_y)) > 0.01) return false;
116 */
117 }
118
119 return true;
120}
121
123{
124 if(cell->getSampling()<21) {
125 m_deta.add(cell->deta());
126 m_dphi.add(cell->dphi());
127 m_mindeta=TMath::Min(cell->deta(),m_mindeta);
128 m_mindphi=TMath::Min(cell->dphi(),m_mindphi);
129
130 m_eta_correction.add(cell->eta_raw()-cell->eta());
131 m_phi_correction.add(cell->phi_raw()-cell->phi());
132 m_mineta_correction=TMath::Min(cell->eta_raw()-cell->eta() , m_mineta_correction);
133 m_maxeta_correction=TMath::Max(cell->eta_raw()-cell->eta() , m_maxeta_correction);
134 m_minphi_correction=TMath::Min(cell->phi_raw()-cell->phi() , m_minphi_correction);
135 m_maxphi_correction=TMath::Max(cell->phi_raw()-cell->phi() , m_maxphi_correction);
136
137 m_mineta=TMath::Min(cell->eta()-cell->deta()/2,m_mineta);
138 m_maxeta=TMath::Max(cell->eta()+cell->deta()/2,m_maxeta);
139 m_minphi=TMath::Min(cell->phi()-cell->dphi()/2,m_minphi);
140 m_maxphi=TMath::Max(cell->phi()+cell->dphi()/2,m_maxphi);
141
142 m_mineta_raw=TMath::Min(cell->eta_raw()-cell->deta()/2,m_mineta_raw);
143 m_maxeta_raw=TMath::Max(cell->eta_raw()+cell->deta()/2,m_maxeta_raw);
144 m_minphi_raw=TMath::Min(cell->phi_raw()-cell->dphi()/2,m_minphi_raw);
145 m_maxphi_raw=TMath::Max(cell->phi_raw()+cell->dphi()/2,m_maxphi_raw);
146 } else {
147 m_deta.add(cell->dx());
148 m_dphi.add(cell->dy());
149 m_mindeta=TMath::Min(cell->dx(),m_mindeta);
150 m_mindphi=TMath::Min(cell->dy(),m_mindphi);
151
152 m_eta_correction.add(cell->x_raw()-cell->x());
153 m_phi_correction.add(cell->y_raw()-cell->y());
154 m_mineta_correction=TMath::Min(cell->x_raw()-cell->x() , m_mineta_correction);
155 m_maxeta_correction=TMath::Max(cell->x_raw()-cell->x() , m_maxeta_correction);
156 m_minphi_correction=TMath::Min(cell->y_raw()-cell->y() , m_minphi_correction);
157 m_maxphi_correction=TMath::Max(cell->y_raw()-cell->y() , m_maxphi_correction);
158
159 m_mineta=TMath::Min(cell->x()-cell->dx()/2,m_mineta);
160 m_maxeta=TMath::Max(cell->x()+cell->dx()/2,m_maxeta);
161 m_minphi=TMath::Min(cell->y()-cell->dy()/2,m_minphi);
162 m_maxphi=TMath::Max(cell->y()+cell->dy()/2,m_maxphi);
163
164 m_mineta_raw=TMath::Min(cell->x_raw()-cell->dx()/2,m_mineta_raw);
165 m_maxeta_raw=TMath::Max(cell->x_raw()+cell->dx()/2,m_maxeta_raw);
166 m_minphi_raw=TMath::Min(cell->y_raw()-cell->dy()/2,m_minphi_raw);
167 m_maxphi_raw=TMath::Max(cell->y_raw()+cell->dy()/2,m_maxphi_raw);
168 }
169 m_cells[cell->identify()]=cell;
170}
171
172bool CaloGeometryLookup::index_range_adjust(int& ieta,int& iphi) const
173{
174 while(iphi<0) {iphi+=m_cell_grid_phi;};
175 while(iphi>=m_cell_grid_phi) {iphi-=m_cell_grid_phi;};
176 if(ieta<0) {
177 ieta=0;
178 return false;
179 }
180 if(ieta>=m_cell_grid_eta) {
181 ieta=m_cell_grid_eta-1;
182 return false;
183 }
184 return true;
185}
186
188{
189 if(size()==0) return;
190 t_cellmap::iterator ic=m_cells.begin();
191 const CaloDetDescrElement* refcell=ic->second;
192 int sampling=refcell->getSampling();
193 if(sampling<21) {
194 double rneta=neta_double()-neta();
195 double rnphi=nphi_double()-nphi();
196 if(TMath::Abs(rneta)>0.05 || TMath::Abs(rnphi)>0.05) {
197 cout<<"ERROR: can't map cells into a regular grid for sampling "<<sampling<<", map index "<<index()<<endl;
198 return;
199 }
204 } else {
205 //double rnx=nx_double()-nx();
206 //double rny=ny_double()-ny();
207 //if(TMath::Abs(rnx)>0.05 || TMath::Abs(rny)>0.05) {
208 // cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": mapping cells into a regular grid, although cells don't fit well"<<endl;
209 //}
210
211 m_cell_grid_eta=TMath::Nint(TMath::Ceil(nx_double()/m_xy_grid_adjustment_factor));
212 m_cell_grid_phi=TMath::Nint(TMath::Ceil(ny_double()/m_xy_grid_adjustment_factor));
215 }
216
218 for(int ieta=0;ieta<m_cell_grid_eta;++ieta) {
219 m_cell_grid[ieta].resize(m_cell_grid_phi,nullptr);
220 }
221
222 for(ic=m_cells.begin();ic!=m_cells.end();++ic) {
223 refcell=ic->second;
224 int ieta,iphi;
225 if(sampling<21) {
226 ieta=raw_eta_position_to_index(refcell->eta_raw());
227 iphi=raw_phi_position_to_index(refcell->phi_raw());
228 } else {
229 ieta=raw_eta_position_to_index(refcell->x_raw());
230 iphi=raw_phi_position_to_index(refcell->y_raw());
231 }
232 if(index_range_adjust(ieta,iphi)) {
233 if(m_cell_grid[ieta][iphi]) {
234 cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": Already cell found at pos ("<<ieta<<","<<iphi<<"): id="
235 <<m_cell_grid[ieta][iphi]->identify()<<"->"<<refcell->identify()<<endl;
236 cout<<" x="<<m_cells[m_cell_grid[ieta][iphi]->identify()]->x_raw()<<" -> "<<refcell->x_raw();
237 cout<<" , y="<<m_cells[m_cell_grid[ieta][iphi]->identify()]->y_raw()<<" -> "<<refcell->y_raw();
238 cout<<" mindx="<<m_deta_double<<" mindy="<<m_dphi_double<<endl;
239 cout<<endl;
240 }
241 m_cell_grid[ieta][iphi]=refcell;
242 } else {
243 cout<<"Grid: Sampling "<<sampling<<"_"<<index()<<": Cell at pos ("<<ieta<<","<<iphi<<"): id="
244 <<refcell->identify()<<" outside eta range"<<endl;
245 }
246 }
247}
248
249float CaloGeometryLookup::calculate_distance_eta_phi(const CaloDetDescrElement* DDE,float eta,float phi,float& dist_eta0,float& dist_phi0) const
250{
251 dist_eta0=(eta - DDE->eta())/m_deta_double;
252 dist_phi0=(TVector2::Phi_mpi_pi(phi - DDE->phi()))/m_dphi_double;
253 float abs_dist_eta0=TMath::Abs(dist_eta0);
254 float abs_dist_phi0=TMath::Abs(dist_phi0);
255 return TMath::Max(abs_dist_eta0,abs_dist_phi0)-0.5;
256}
257
258const CaloDetDescrElement* CaloGeometryLookup::getDDE(float eta,float phi,float* distance,int* steps)
259{
260 float dist = 0;
261 const CaloDetDescrElement* bestDDE=nullptr;
262 if(!distance) distance=&dist;
263 *distance=+10000000;
264 int intsteps=0;
265 if(!steps) steps=&intsteps;
266
267 float best_eta_corr=m_eta_correction;
268 float best_phi_corr=m_phi_correction;
269
270 float raw_eta=eta+best_eta_corr;
271 float raw_phi=phi+best_phi_corr;
272
273 int ieta=raw_eta_position_to_index(raw_eta);
274 int iphi=raw_phi_position_to_index(raw_phi);
275 index_range_adjust(ieta,iphi);
276
277 const CaloDetDescrElement* newDDE=m_cell_grid[ieta][iphi];
278 float bestdist=+10000000;
279 ++(*steps);
280 int nsearch=0;
281 while(newDDE && nsearch<3) {
282 float dist_eta0,dist_phi0;
283 *distance=calculate_distance_eta_phi(newDDE,eta,phi,dist_eta0,dist_phi0);
284 bestDDE=newDDE;
285 bestdist=*distance;
286
288 cout<<"CaloGeometryLookup::getDDE, search iteration "<<nsearch<<endl;
289 cout<<"cell id="<<newDDE->identify()<<" steps "<<*steps<<" steps, eta="<<eta<<" phi="<<phi<<" dist="<<*distance<<" cell_grid_eta="<<cell_grid_eta()<<" cell_grid_phi="<<cell_grid_phi()<<" m_deta_double="<<m_deta_double<<" m_dphi_double="<<m_dphi_double<<" 1st_eta_corr="<<best_eta_corr<<" 1st_phi_corr="<<best_phi_corr<<endl;
290 cout<<" input sampling="<<newDDE->getSampling()<<" eta="<<newDDE->eta()<<" eta_raw="<<newDDE->eta_raw()<<" deta="<<newDDE->deta()<<" ("<<(newDDE->eta_raw()-newDDE->eta())/newDDE->deta()<<") phi="<<newDDE->phi()<<" phi_raw="<<newDDE->phi_raw()<<" dphi="<<newDDE->dphi()<<" ("<<(newDDE->phi_raw()-newDDE->phi())/newDDE->dphi()<<")"<<endl;
291 cout<<" ieta="<<ieta<<" iphi="<<iphi<<endl;
292 cout<<" dist_eta0="<<dist_eta0<<" ,dist_phi0="<<dist_phi0<<endl;
293 }
294
295 if(*distance<0) return newDDE;
296 //correct ieta and iphi by the observed difference to the hit cell
297 ieta+=TMath::Nint(dist_eta0);
298 iphi+=TMath::Nint(dist_phi0);
299 index_range_adjust(ieta,iphi);
300 const CaloDetDescrElement* oldDDE=newDDE;
301 newDDE=m_cell_grid[ieta][iphi];
302 ++(*steps);
303 ++nsearch;
304 if(oldDDE==newDDE) break;
305 }
306 if(CaloGeometry::m_debug && !newDDE) {
307 cout<<"CaloGeometryLookup::getDDE, direct search ieta="<<ieta<<" iphi="<<iphi<<" is empty"<<endl;
308 }
309 float minieta=ieta+TMath::Nint(TMath::Floor(m_mineta_correction/cell_grid_eta()));
310 float maxieta=ieta+TMath::Nint(TMath::Ceil(m_maxeta_correction/cell_grid_eta()));
311 float miniphi=iphi+TMath::Nint(TMath::Floor(m_minphi_correction/cell_grid_phi()));
312 float maxiphi=iphi+TMath::Nint(TMath::Ceil(m_maxphi_correction/cell_grid_phi()));
313 if(minieta<0) minieta=0;
314 if(maxieta>=m_cell_grid_eta) maxieta=m_cell_grid_eta-1;
315 for(int iieta=minieta;iieta<=maxieta;++iieta) {
316 for(int iiphi=miniphi;iiphi<=maxiphi;++iiphi) {
317 ieta=iieta;
318 iphi=iiphi;
319 index_range_adjust(ieta,iphi);
320 newDDE=m_cell_grid[ieta][iphi];
321 ++(*steps);
322 if(newDDE) {
323 float dist_eta0,dist_phi0;
324 *distance=calculate_distance_eta_phi(newDDE,eta,phi,dist_eta0,dist_phi0);
325
327 cout<<"CaloGeometryLookup::getDDE, windows search! minieta="<<m_mineta_correction/cell_grid_eta()<<" maxieta="<<m_maxeta_correction/cell_grid_eta()<<" miniphi="<<m_minphi_correction/cell_grid_phi()<<" maxiphi="<<m_maxphi_correction/cell_grid_phi()<<endl;
328 cout<<"cell id="<<newDDE->identify()<<" steps "<<*steps<<" steps, eta="<<eta<<" phi="<<phi<<" dist="<<*distance<<endl;
329 cout<<" input sampling="<<newDDE->getSampling()<<" eta="<<newDDE->eta()<<" eta_raw="<<newDDE->eta_raw()<<" deta="<<newDDE->deta()<<" ("<<(newDDE->eta_raw()-newDDE->eta())/newDDE->deta()<<") phi="<<newDDE->phi()<<" phi_raw="<<newDDE->phi_raw()<<" dphi="<<newDDE->dphi()<<" ("<<(newDDE->phi_raw()-newDDE->phi())/newDDE->dphi()<<")"<<endl;
330 cout<<" ieta="<<ieta<<" iphi="<<iphi<<endl;
331 cout<<" dist_eta0="<<dist_eta0<<" ,dist_phi0="<<dist_phi0<<endl;
332 }
333
334 if(*distance<0) return newDDE;
335 if(*distance<bestdist) {
336 bestDDE=newDDE;
337 bestdist=*distance;
338 }
339 } else if(CaloGeometry::m_debug) {
340 cout<<"CaloGeometryLookup::getDDE, windows search ieta="<<ieta<<" iphi="<<iphi<<" is empty"<<endl;
341 }
342 }
343 }
344 *distance=bestdist;
345 return bestDDE;
346}
347
348
349
350/*
351void CaloGeometryLookup::CalculateTransformation()
352{
353 gROOT->cd();
354
355 TTree* tmap = new TTree( "mapping" , "mapping" );
356
357 float eta,phi,Deta_raw,Dphi_raw;
358 tmap->Branch("eta", &eta,"eta/F");
359 tmap->Branch("phi", &phi,"phi/F");
360 tmap->Branch("Deta_raw", &Deta_raw,"Deta_raw/F");
361 tmap->Branch("Dphi_raw", &Dphi_raw,"Dphi_raw/F");
362
363 if(m_cells.size()==0) return;
364
365 int sampling=0;
366 for(t_cellmap::iterator ic=m_cells.begin();ic!=m_cells.end();++ic) {
367 CaloGeoDetDescrElement* refcell=ic->second;
368 sampling=refcell->getSampling();
369 if(sampling<21) {
370 eta=refcell->eta();
371 phi=refcell->phi();
372 Deta_raw=refcell->eta_raw()-eta;
373 Dphi_raw=refcell->phi_raw()-phi;
374 } else {
375 eta=refcell->x();
376 phi=refcell->y();
377 Deta_raw=refcell->x_raw()-eta;
378 Dphi_raw=refcell->y_raw()-phi;
379 }
380 tmap->Fill();
381 tmap->Fill(); //Fill twice to have all events and training and test tree
382 }
383 tmap->Print();
384
385 TString outfileName( Form("Mapping%d_%d.root",sampling,m_index) );
386 TFile* outputFile = new TFile( outfileName, "RECREATE" );
387 //TFile* outputFile = new TMemFile( outfileName, "RECREATE" );
388
389 TMVA::Factory *factory = new TMVA::Factory( Form("Mapping%d_%d.root",sampling,m_index) , outputFile, "!V:!Silent:Color:DrawProgressBar" );
390
391 factory->AddVariable( "eta", "calo eta", "1", 'F' );
392 factory->AddVariable( "phi", "calo phi", "1", 'F' );
393 factory->AddTarget( "Deta_raw" );
394 factory->AddTarget( "Dphi_raw" );
395
396 factory->AddRegressionTree( tmap );
397
398 //factory->PrepareTrainingAndTestTree( "" , Form("nTrain_Regression=%d:nTest_Regression=%d:SplitMode=Random:NormMode=NumEvents:!V",(int)m_cells.size(),(int)m_cells.size()) );
399 factory->PrepareTrainingAndTestTree( "" , "nTrain_Regression=0:nTest_Regression=0:SplitMode=Alternate:NormMode=NumEvents:!V" );
400
401 factory->BookMethod( TMVA::Types::kMLP, "MLP", "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N+5:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-6:ConvergenceTests=15:!UseRegulator" );
402
403 cout<<"=== Start trainging ==="<<endl;
404 // Train MVAs using the set of training events
405 factory->TrainAllMethods();
406
407 cout<<"=== Start testing ==="<<endl;
408 // ---- Evaluate all MVAs using the set of test events
409 factory->TestAllMethods();
410
411 cout<<"=== Start evaluation ==="<<endl;
412 // ----- Evaluate and compare performance of all configured MVAs
413 factory->EvaluateAllMethods();
414
415 outputFile->Close();
416
417 delete factory;
418
419 delete tmap;
420}
421*/
const boost::regex ref(r_ef)
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
This class groups all DetDescr information related to a CaloCell.
CaloCell_ID::CaloSample getSampling() const
cell sampling
Identifier identify() const override final
cell identifier
const MeanAndRMS & deta()
std::vector< std::vector< const CaloDetDescrElement * > > m_cell_grid
const MeanAndRMS & dphi()
bool IsCompatible(const CaloDetDescrElement *cell)
int raw_phi_position_to_index(float phi_raw) const
t_cellmap::size_type size() const
void add(const CaloDetDescrElement *cell)
virtual const CaloDetDescrElement * getDDE(float eta, float phi, float *distance=0, int *steps=0)
bool index_range_adjust(int &ieta, int &iphi) const
bool has_overlap(CaloGeometryLookup *ref)
int raw_eta_position_to_index(float eta_raw) const
float calculate_distance_eta_phi(const CaloDetDescrElement *DDE, float eta, float phi, float &dist_eta0, float &dist_phi0) const
virtual ~CaloGeometryLookup()
void merge_into_ref(CaloGeometryLookup *ref)
static std::atomic< bool > m_debug
static const Identifier m_debug_identify
STL namespace.