]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliJetGrid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Process codes corrected/completed.
[u/mrichter/AliRoot.git] / JETAN / AliJetGrid.cxx
CommitLineData
4a01bb2c 1/**************************************************************************
2 * Copyright(c) 2001-2002, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: AliJetGrid.cxx,v 1.0 07/09/2006 */
17
18//=========================================================================
19// *** 07 September 2006
20// This class allows to fill a really general grid in (eta,phi)
21// Two types of grid can be setted :
22// if fGrid==0 -> the complete acceptance of a rectangular detector acceptance is filled
23// if fGrid==1 -> the previous grid - an other rectangular detector acceptance is filled
24// A (eta,phi(rad)) position can be extracted from an index value without looping over all entries
25// An index position can be extracted from a (eta,phi(rad)) position without looping over all entries
26//
27// How to run it :
28// > AliJetGrid *grid = new AliJetGrid((NbinPhi-1),(NbinEta-1),PhiMin,PhiMax,EtaMin,EtaMax);
29// > grid->SetGridType(...); // 0 or 1 for the moment
30// > grid->SetMatrixIndexes();
31// > grid->SetIndexIJ();
32//
33// Author : magali.estienne@ires.in2p3.fr
34//=========================================================================
35
36// Standard headers
37#include <Riostream.h>
38// Root headers
37b09b91 39#include <TMath.h>
4a01bb2c 40#include <TMatrixD.h>
41#include <TArrayD.h>
42#include <TArrayI.h>
43// AliRoot headers
44#include "AliJetGrid.h"
45
87f40625 46ClassImp(AliJetGrid)
4a01bb2c 47
4a01bb2c 48
b92e2ccf 49//__________________________________________________________
50AliJetGrid::AliJetGrid():
51 fGrid(0),
52 fNphi(0),
53 fNeta(0),
54 fPhi(0),
55 fEta(0),
56 fIndex(0),
57 fIndexI(0),
58 fIndexJ(0),
59 fPhiMin(0),
60 fPhiMax(0),
61 fEtaMin(0),
62 fEtaMax(0),
63 fEtaBinInTPCAcc(0),
64 fPhiBinInTPCAcc(0),
65 fEtaBinInEMCalAcc(0),
66 fPhiBinInEMCalAcc(0),
67 fNbinEta(0),
68 fNbinPhi(0),
69 fMaxPhi(0),
70 fMinPhi(0),
71 fMaxEta(0),
72 fMinEta(0),
73 fDebug(1)
74 {
4a01bb2c 75 // Default constructor
4a01bb2c 76}
77
78//__________________________________________________________
b92e2ccf 79AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax):
80 fGrid(0),
81 fNphi(nphi),
82 fNeta(neta),
83 fPhi(0),
84 fEta(0),
85 fIndex(0),
86 fIndexI(0),
87 fIndexJ(0),
88 fPhiMin(0),
89 fPhiMax(0),
90 fEtaMin(0),
91 fEtaMax(0),
92 fEtaBinInTPCAcc(0),
93 fPhiBinInTPCAcc(0),
94 fEtaBinInEMCalAcc(0),
95 fPhiBinInEMCalAcc(0),
96 fNbinEta(0),
97 fNbinPhi(0),
98 fMaxPhi(phiMax),
99 fMinPhi(phiMin),
100 fMaxEta(etaMax),
101 fMinEta(etaMin),
102 fDebug(1)
103 {
4a01bb2c 104 // Standard constructor
b92e2ccf 105 fPhi = new TArrayD(fNphi+1);
106 fEta = new TArrayD(fNeta+1);
4a01bb2c 107 fIndexI = new TArrayI((fNeta+1)*(fNphi+1)+1);
108 fIndexJ = new TArrayI((fNeta+1)*(fNphi+1)+1);
b92e2ccf 109
4a01bb2c 110 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = (phiMax-phiMin)/fNphi*i+phiMin;
111 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = (etaMax-etaMin)/fNeta*i+etaMin;
b92e2ccf 112
4a01bb2c 113 if(fDebug > 3){
b92e2ccf 114 for(Int_t i=0; i<(fNphi+1); i++) cout << (*fPhi)[i] << endl;
115 for(Int_t i=0; i<(fNeta+1); i++) cout << (*fEta)[i] << endl;
4a01bb2c 116 }
b92e2ccf 117
4a01bb2c 118 fIndex = new TMatrixD(fNphi+1,fNeta+1);
b92e2ccf 119
4a01bb2c 120}
121
122//__________________________________________________________
123AliJetGrid::AliJetGrid(const AliJetGrid& grid):TNamed(grid) {
124
125 // Copy constructor
126
127 fNphi = grid.fNphi;
128 fNeta = grid.fNeta;
129 fPhiMin = grid.fPhiMin;
130 fPhiMax = grid.fPhiMax;
131 fEtaMin = grid.fEtaMin;
132 fEtaMax = grid.fEtaMax;
133 fEtaBinInTPCAcc = grid.fEtaBinInTPCAcc;
134 fPhiBinInTPCAcc = grid.fPhiBinInTPCAcc;
135 fEtaBinInEMCalAcc = grid.fEtaBinInEMCalAcc;
136 fPhiBinInEMCalAcc = grid.fPhiBinInEMCalAcc;
137 fNbinPhi = grid.fNbinPhi;
138 fMaxPhi = grid.fMaxPhi;
139 fMinPhi = grid.fMinPhi;
140 fMaxEta = grid.fMaxEta;
141 fMinEta = grid.fMinEta;
142
143 fPhi = new TArrayD(fNphi+1);
144 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = grid.fPhi->At(i);
145 fEta = new TArrayD(fNeta+1);
146 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = grid.fEta->At(i);
147
148 fIndex = new TMatrixD(fNphi+1,fNeta+1);
149 for(Int_t i=0; i<fNphi+1; i++) {
150 for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*grid.fIndex)(i,j);
151 }
152}
153
154//__________________________________________________________
155AliJetGrid::~AliJetGrid() {
156
157 // Destructor
158 delete fPhi;
159 delete fEta;
160 delete fIndexI;
161 delete fIndexJ;
162 delete fIndex;
163}
164
165//__________________________________________________________
166void AliJetGrid::InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal)
d980dfdb 167{
168// To set initial parameters
4a01bb2c 169
170 fPhiMin = phiMinCal; // rad
171 fPhiMax = phiMaxCal; // rad
172 fEtaMin = etaMinCal;
173 fEtaMax = etaMaxCal;
174 fNbinPhi = static_cast<int>(fPhiMin/(2*TMath::Pi()/fNphi));
175
176 // Define some binning numbers
177 if(fGrid==0){
178 for(Int_t i=0; i<fNphi+1; i++) fPhiBinInTPCAcc++;
179 fPhiBinInEMCalAcc = 0;
180
181 for(Int_t i=0; i<fNeta+1; i++) fEtaBinInTPCAcc++;
182 fEtaBinInEMCalAcc = 0;
183 }
184
185 if(fGrid==1){
186 for(Int_t i=0; i<fNphi+1; i++) {
187 fPhiBinInTPCAcc++;
188 if(fPhi->At(i) >= fPhiMin &&
189 fPhi->At(i) <= fPhiMax)
190 fPhiBinInEMCalAcc++;
191 }
192 for(Int_t i=0; i<fNeta+1; i++) {
193 fEtaBinInTPCAcc++;
194 if((fEta->At(i) >= fEtaMin) &&
195 (fEta->At(i) <= fEtaMax))
196 fEtaBinInEMCalAcc++;
197 }
198 }
199
200}
201
202//__________________________________________________________
203TArrayD* AliJetGrid::GetArrayEta()
d980dfdb 204{
205// Returns an array with the eta points
4a01bb2c 206
207 return fEta;
208
209}
210
211//__________________________________________________________
212TArrayD* AliJetGrid::GetArrayPhi()
d980dfdb 213{
214// Returns an array with the phi points
4a01bb2c 215
216 return fPhi;
217
218}
219
220//__________________________________________________________
221TMatrixD* AliJetGrid::GetIndexObject()
d980dfdb 222{
223// Returns a pointer to the matrix
4a01bb2c 224
225 return fIndex;
226
227}
228
229//__________________________________________________________
230void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi,
d980dfdb 231 Float_t &mineta, Float_t &maxeta) const
232{
233// Returns EMCAL acceptance initially setted
4a01bb2c 234
235 nphi = fNphi;
236 neta = fNeta;
237 minphi = fPhiMin;
238 maxphi = fPhiMax;
239 mineta = fEtaMin;
240 maxeta = fEtaMax;
241}
242
243//__________________________________________________________
244void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc,
d980dfdb 245 Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi) const
246{
247// Returns number of bins in TPC and EMCAL geometry
4a01bb2c 248
249 etabintpc = fEtaBinInTPCAcc;
250 phibintpc = fPhiBinInTPCAcc;
251 etabinemc = fEtaBinInEMCalAcc;
252 phibinemc = fPhiBinInEMCalAcc;
253 nbinphi = fNbinPhi;
254}
255
256//__________________________________________________________
257Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const
d980dfdb 258{
259// Tells the index value of a corresponding (eta,phi) real position
260// Loop over all entries -> takes time.
261// Used one time at the begining to fill the grids
4a01bb2c 262
263 /* this is how bins are numbered
264
265 in all TPC or in TPC - EMCAL
266 ... ... ... ..
267 ---------------
268 ... ... . 10 | | | 11
269 ---+---+--- ---------------
270 ^ 6 | 7 | 8 or 8 | | | 9
271 | ---+---+--- ---------------
272 3 | 4 | 5 4 | 5 | 6 | 7
273 phi ---+---+--- ---------------
274 0 | 1 | 2 0 | 1 | 2 | 3
275
276
277 eta ->
278 */
279
280 Int_t etaBin=0,phiBin=0,absID=0;
281 Int_t etaBin2=0,etaBin3=0;
282
283 // Fill all the grid in eta/phi (all TPC acceptance)
284 //-----------------------------------------------------
285 if(fGrid==0){
286 if(eta <= fEta->At(0)) {
287 etaBin = 0;
288 } else if(eta >= fEta->At(fNeta)) {
289 etaBin = fNeta;
290 } else {
291 for(Int_t i=0; i<fNeta+1; i++) {
292 if(eta < fEta->At(i)) {
293 etaBin = i-1;
294 break;
295 }
296 }
297 }
298 if(phi <= fPhi->At(0)) {
299 phiBin = 0;
300 } else if(phi >= fPhi->At(fNphi)) {
301 phiBin = fNphi;
302 } else {
303 for(Int_t i=0; i<fNphi+1; i++) {
304 if(phi < fPhi->At(i)) {
305 phiBin = i-1;
306 break;
307 }
308 }
309 }
310
311 // Calculate absolute id
312 absID = phiBin*(fNeta+1) + etaBin;
313
314 }
315
316 // Fill the grid but do not count id in EMCal acceptance
317 //------------------------------------------------------
318 if((eta >= fEtaMin && eta <= fEtaMax) &&
319 (phi >= fPhiMin && phi <= fPhiMax)){
320 etaBin = etaBin2 = etaBin3 = -1;
321 phiBin = -1;
322 }
323
324 if(fGrid == 1){
325 if(phi<fPhiMin) {
326 if(eta <= fEta->At(0)) {
327 etaBin = 0;
328 } else if(eta >= fEta->At(fNeta)) {
329 etaBin = fNeta;
330 } else {
331 for(Int_t i=0; i<fNeta+1; i++) {
332 if(eta < fEta->At(i)) {
333 etaBin = i-1;
334 break;
335 }
336 }
337 }
338 }
339 else {
340 if((phi>=fPhiMin) && (phi<=fPhiMax)) {
341 if(eta <= fEta->At(0)) {
342 etaBin2 = 0;
343 } else if(eta >= fEta->At(fNeta)) {
344 etaBin2 = fNeta-fEtaBinInEMCalAcc;
345 } else {
346 for(Int_t i=0; i<fNeta+1; i++) {
347 if(eta < fEta->At(i) && ((fEta->At(0)<eta && eta<fEtaMin) ||
348 (fEtaMax<eta && eta<fEta->At(fNeta)))) {
349 // cout << "i : " << i << endl;
350 if(eta<fEtaMin)
351 etaBin2 = i-1;
352 else etaBin2 = i-1-fEtaBinInEMCalAcc;
353 break;
354 }
355 }
356 }
357 }
358 else {
359 if(eta <= fEta->At(0)) {
360 etaBin3 = 0;
361 } else if(eta >= fEta->At(fNeta)) {
362 etaBin3 = fNeta;
363 } else {
364 for(Int_t i=0; i<fNeta+1; i++) {
365 if(eta < fEta->At(i)) {
366 etaBin3 = i-1;
367 break;
368 }
369 }
370 }
371 }
372 }
373
374 if(phi <= fPhi->At(0)) {
375 phiBin = 0;
376 } else if(phi >= fPhi->At(fNphi)) {
377 phiBin = fNphi;
378 } else {
379 for(Int_t i=0; i<fNphi+1; i++) {
380 if(phi < fPhi->At(i)) {
381 phiBin = i-1;
382 break;
383 }
384 }
385 }
386
387 // Calculate absID
388 //-------------------------------------------------------
389 if(phi<fPhiMin)
390 absID = phiBin*(fNeta+1) + etaBin;
391 if(phi>=fPhiMin && phi<=fPhiMax) {
392 if(eta >= fEtaMin && eta <= fEtaMax) absID = -1;
393 else{
394 absID = (fNbinPhi+1)*(fNeta+1)
395 + (phiBin-fNbinPhi-1)*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
396 + etaBin2;
397 }
398 }
399 if(phi>fPhiMax)
400 absID = (fNbinPhi+1)*(fNeta+1)
401 + fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
402 + (phiBin-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)
403 + etaBin3;
404
405 } // END OPTION==1
406
407 return absID;
408
409}
410
411//__________________________________________________________
412void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
d980dfdb 413{
414// Get (eta,phi) position for a given index BUT loop over all entries (takes time)
4a01bb2c 415
416 for(Int_t j=0; j<fNphi+1; j++) {
417 for(Int_t i=0; i<fNeta+1; i++) {
418
419 // TPC grid only
420 //-------------------------------------
421 if(fGrid==0) {
422 if(j*(fNeta+1)+i == index) {
423 eta = fEta->At(i);
424 phi = fPhi->At(j);
425 }
426 }
427
428 // TPC-EMCAL grid
429 //-------------------------------------
430 Int_t ii = 0;
431 if(i==0) ii = 0;
432 if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
433 if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
434
435 if(fGrid==1) {
436 if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
437 eta = fEta->At(i);
438 phi = fPhi->At(j);
439 }
440
441 if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
442 ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
443 if(ii==0) {Int_t ind = 0; eta = fEta->At(ind);}
444 else eta = fEta->At(i);
445 phi = fPhi->At(j);
446 }
447
448 if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) &&
449 ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
450 +(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
451 eta = fEta->At(i);
452 phi = fPhi->At(j);
453 }
454 }
455 }
456 }
457}
458
459//__________________________________________________________
460Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta)
d980dfdb 461{
462// Get index value for a (eta,phi) position - Direct value
4a01bb2c 463
464 Int_t ieta = GetIndexJFromEta(eta);
465 Int_t iphi = GetIndexIFromPhi(phi);
466
467 Int_t index = GetMatrixIndex(iphi,ieta);
468
469 if(fDebug>10){
470 cout << "(phi,eta) : " << phi << ", " << eta << endl;
471 cout << "index : " << index << endl;
472 }
473 return index;
474
475}
476
477//__________________________________________________________
478Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
d980dfdb 479{
480// Get eta id
481// Eta discretized
4a01bb2c 482
483 Int_t idEta =0;
484 Double_t temp = (eta+fMaxEta)/(fMaxEta-fMinEta)*fNeta;
485 if(fDebug>20)
486 {
487 cout << "eta : " << eta << endl;
488 cout << "fMaxEta : " << fMaxEta << endl;
489 cout << "fMinEta : " << fMinEta << endl;
490 cout << "fNeta : " << fNeta << endl;
491 cout << "temp eta before cast : " << temp << endl;
492 }
493 idEta = static_cast<Int_t>(temp+0.5);
494 if(fDebug>20) cout << "temp eta after cast : " << idEta << endl;
495 return idEta;
496}
497//__________________________________________________________
498Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
d980dfdb 499{
500// Get phi id
501// Phi discretized
4a01bb2c 502
503 Int_t idPhi = 0;
504 Double_t temp = 0.;
505 if(fMinPhi==0) temp = phi/(fMaxPhi-fMinPhi)*fNphi;
506 else temp = (phi-fMinPhi)/(fMaxPhi-fMinPhi)*fNphi;
507
508 if(fDebug>20)
509 {
510 cout << "phi : " << phi << endl;
511 cout << "fMaxPhi : " << fMaxPhi << endl;
512 cout << "fMinPhi : " << fMinPhi << endl;
513 cout << "fNphi : " << fNphi << endl;
514 cout << "temp phi before cast : " << temp << endl;
515 }
516 idPhi = static_cast<Int_t>(temp+0.5);
517 if(fDebug>20) cout << "temp phi after cast : " << idPhi << endl;
518 return idPhi;
519
520
521}
522
523//__________________________________________________________
524void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par)
d980dfdb 525{
526// Allows to set parameters using only one index (if fGrid==0) !!
527// Not used !
4a01bb2c 528
529 Int_t iphi = (Int_t)i/fNeta;
530 Int_t ieta = i-iphi*fNeta;
531 SetMatrixIndex(iphi,ieta,par);
532
533 return;
534}
535
536//__________________________________________________________
537void AliJetGrid::SetMatrixIndexes()
d980dfdb 538{
539// Fill the final matrix object with the corresponding index in eta/phi
4a01bb2c 540
541 for(Int_t i=0; i<fNphi+1; i++){
542 for(Int_t j=0; j<fNeta+1; j++){
543 (*fIndex)(i,j) = GetIndexFromEtaPhi(fPhi->At(i),fEta->At(j))+1;
544 if(fDebug>2){
545 cout << "(*fIndex)(" << i << "," << j << ") : " << (*fIndex)(i,j) <<
546 ", phi : " << fPhi->At(i) << ", eta : " << fEta->At(j) << endl;
547 }
548 }
549 }
550 printf("##########################################################\n");
551 printf("TMatrix object filled !\n");
552 printf("Size of the object : phi x eta = (fNphi+1) x (fNeta+1) = %d\n",(fNphi+1)*(fNeta+1));
553 printf("##########################################################\n");
554}
555
556//__________________________________________________________
557void AliJetGrid::SetIndexIJ()
d980dfdb 558{
559// Set the grid index
4a01bb2c 560
561 for(Int_t i=0; i<fNphi+1; i++){
562 for(Int_t j=0; j<fNeta+1; j++){
563 Int_t id = static_cast<Int_t>((*fIndex)(i,j));
564
565 if(id!=-1)
566 {
567 (*fIndexI)[id] = i;
568 (*fIndexJ)[id] = j;
569 }
570 }
571 }
572
573 printf("##########################################################\n");
574 printf(" In SetIndexIJ - Grid indexes setted !\n");
575 printf("##########################################################\n");
576}
577
578//__________________________________________________________
d980dfdb 579void AliJetGrid::GetIJFromIndex(Int_t index, Int_t i, Int_t j) const
580{
581// Returns i position id of eta and j position id of phi for a given grid index
4a01bb2c 582 i = (*fIndexI)[index];
583 j = (*fIndexJ)[index];
584}
585
586//__________________________________________________________
587void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
d980dfdb 588{
589// Returns eta, phi values for a given grid index
4a01bb2c 590
591 phi = fPhi->At((*fIndexI)[index]);
592 eta = fEta->At((*fIndexJ)[index]);
593}
594
595//__________________________________________________________
596Int_t AliJetGrid::GetNEntries()
d980dfdb 597{
598// Returns the number of entries of the grid
4a01bb2c 599
600 if(fDebug>20){
601 cout << "fMaxPhi : " << fMaxPhi << endl;
602 cout << "fMaxEta : " << fMaxEta << endl;
603 }
604
605 Int_t indexNum = GetIndex(fMaxPhi,fMaxEta);
606 if(fDebug>20) cout << "indexNum : " << indexNum << endl;
607 return indexNum;
608
609}
610
611//__________________________________________________________
612Int_t AliJetGrid::GetNEntries2()
d980dfdb 613{
614// Returns the number of entries of the grid
4a01bb2c 615
616 Int_t indexNum = GetIndex(fMaxPhi-1.,fMaxEta-0.5);
617 if(fDebug>20) cout << "indexNum : " << indexNum << endl;
618 return indexNum;
619
620}
621
622
623
624
625
626