1 /**************************************************************************
2 * Copyright(c) 2001-2002, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliJetGrid.cxx,v 1.0 07/09/2006 */
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
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();
33 // Author : magali.estienne@subatech.in2p3.fr
34 //=========================================================================
37 #include <Riostream.h>
44 #include "AliJetGrid.h"
48 //__________________________________________________________
49 AliJetGrid::AliJetGrid():
74 // Default constructor
77 //__________________________________________________________
78 AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax):
103 // Standard constructor
104 fPhi = new TArrayD(fNphi+1);
105 fEta = new TArrayD(fNeta+1);
106 fIndexI = new TArrayI((fNeta+1)*(fNphi+1)+1);
107 fIndexJ = new TArrayI((fNeta+1)*(fNphi+1)+1);
109 for(Int_t i=0; i<fNphi+1; i++) {
110 if(fNphi!=0) (*fPhi)[i] = (phiMax-phiMin)/fNphi*i+phiMin;
111 else (*fPhi)[i] = phiMin+(phiMax-phiMin)/2;
113 for(Int_t i=0; i<fNeta+1; i++) {
114 if(fNeta!=0) (*fEta)[i] = (etaMax-etaMin)/fNeta*i+etaMin;
115 else (*fEta)[i] = etaMin+(etaMax-etaMin)/2;
119 for(Int_t i=0; i<(fNphi+1); i++) cout << (*fPhi)[i] << endl;
120 for(Int_t i=0; i<(fNeta+1); i++) cout << (*fEta)[i] << endl;
123 fIndex = new TMatrixD(fNphi+1,fNeta+1);
127 //__________________________________________________________
128 AliJetGrid::AliJetGrid(const AliJetGrid& grid) :
136 fIndexI(grid.fIndexI),
137 fIndexJ(grid.fIndexJ),
138 fPhiMin(grid.fPhiMin),
139 fPhiMax(grid.fPhiMax),
140 fEtaMin(grid.fEtaMin),
141 fEtaMax(grid.fEtaMax),
142 fEtaBinInTPCAcc(grid.fEtaBinInTPCAcc),
143 fPhiBinInTPCAcc(grid.fPhiBinInTPCAcc),
144 fEtaBinInEMCalAcc(grid.fEtaBinInEMCalAcc),
145 fPhiBinInEMCalAcc(grid.fPhiBinInEMCalAcc),
146 fNbinEta(grid.fNbinEta),
147 fNbinPhi(grid.fNbinPhi),
148 fMaxPhi(grid.fMaxPhi),
149 fMinPhi(grid.fMinPhi),
150 fMaxEta(grid.fMaxEta),
151 fMinEta(grid.fMinEta),
157 fPhi = new TArrayD(fNphi+1);
158 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = grid.fPhi->At(i);
159 fEta = new TArrayD(fNeta+1);
160 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = grid.fEta->At(i);
162 fIndex = new TMatrixD(fNphi+1,fNeta+1);
163 for(Int_t i=0; i<fNphi+1; i++) {
164 for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*grid.fIndex)(i,j);
169 AliJetGrid& AliJetGrid::operator=(const AliJetGrid& other)
178 fIndexI = other.fIndexI;
179 fIndexJ = other.fIndexJ;
180 fPhiMin = other.fPhiMin;
181 fPhiMax = other.fPhiMax;
182 fEtaMin = other.fEtaMin;
183 fEtaMax = other.fEtaMax;
184 fEtaBinInTPCAcc = other.fEtaBinInTPCAcc;
185 fPhiBinInTPCAcc = other.fPhiBinInTPCAcc;
186 fEtaBinInEMCalAcc = other.fEtaBinInEMCalAcc;
187 fPhiBinInEMCalAcc = other.fPhiBinInEMCalAcc;
188 fNbinEta = other.fNbinEta;
189 fNbinPhi = other.fNbinPhi;
190 fMaxPhi = other.fMaxPhi;
191 fMinPhi = other.fMinPhi;
192 fMaxEta = other.fMaxEta;
193 fMinEta = other.fMinEta;
194 fDebug = other.fDebug;
195 fPhi = new TArrayD(fNphi+1);
196 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = other.fPhi->At(i);
197 fEta = new TArrayD(fNeta+1);
198 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = other.fEta->At(i);
200 fIndex = new TMatrixD(fNphi+1,fNeta+1);
201 for(Int_t i=0; i<fNphi+1; i++) {
202 for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*other.fIndex)(i,j);
207 //__________________________________________________________
208 AliJetGrid::~AliJetGrid() {
218 //__________________________________________________________
219 void AliJetGrid::InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal)
221 // To set initial parameters
223 fPhiMin = phiMinCal; // rad
224 fPhiMax = phiMaxCal; // rad
227 fNbinPhi = static_cast<int>(fPhiMin/(2*TMath::Pi()/fNphi));
229 // Define some binning numbers
231 for(Int_t i=0; i<fNphi+1; i++) fPhiBinInTPCAcc++;
232 fPhiBinInEMCalAcc = 0;
234 for(Int_t i=0; i<fNeta+1; i++) fEtaBinInTPCAcc++;
235 fEtaBinInEMCalAcc = 0;
239 for(Int_t i=0; i<fNphi+1; i++) {
241 if(fPhi->At(i) >= fPhiMin &&
242 fPhi->At(i) <= fPhiMax)
245 for(Int_t i=0; i<fNeta+1; i++) {
247 if((fEta->At(i) >= fEtaMin) &&
248 (fEta->At(i) <= fEtaMax))
255 //__________________________________________________________
256 TArrayD* AliJetGrid::GetArrayEta()
258 // Returns an array with the eta points
264 //__________________________________________________________
265 TArrayD* AliJetGrid::GetArrayPhi()
267 // Returns an array with the phi points
273 //__________________________________________________________
274 TMatrixD* AliJetGrid::GetIndexObject()
276 // Returns a pointer to the matrix
282 //__________________________________________________________
283 void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi,
284 Float_t &mineta, Float_t &maxeta) const
286 // Returns EMCAL acceptance initially setted
296 //__________________________________________________________
297 void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc,
298 Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi) const
300 // Returns number of bins in TPC and EMCAL geometry
302 etabintpc = fEtaBinInTPCAcc;
303 phibintpc = fPhiBinInTPCAcc;
304 etabinemc = fEtaBinInEMCalAcc;
305 phibinemc = fPhiBinInEMCalAcc;
309 //__________________________________________________________
310 Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const
312 // Tells the index value of a corresponding (eta,phi) real position
313 // Loop over all entries -> takes time.
314 // Used one time at the begining to fill the grids
316 /* this is how bins are numbered
318 in all TPC or in TPC - EMCAL
321 ... ... . 10 | | | 11
322 ---+---+--- ---------------
323 ^ 6 | 7 | 8 or 8 | | | 9
324 | ---+---+--- ---------------
325 3 | 4 | 5 4 | 5 | 6 | 7
326 phi ---+---+--- ---------------
327 0 | 1 | 2 0 | 1 | 2 | 3
333 Int_t etaBin=0,phiBin=0,absID=0;
334 Int_t etaBin2=0,etaBin3=0;
336 // Fill all the grid in eta/phi (all TPC acceptance)
337 //-----------------------------------------------------
339 if(eta <= fEta->At(0)) {
341 } else if(eta >= fEta->At(fNeta)) {
344 for(Int_t i=0; i<fNeta+1; i++) {
345 if(eta < fEta->At(i)) {
351 if(phi <= fPhi->At(0)) {
353 } else if(phi >= fPhi->At(fNphi)) {
356 for(Int_t i=0; i<fNphi+1; i++) {
357 if(phi < fPhi->At(i)) {
364 // Calculate absolute id
365 absID = phiBin*(fNeta+1) + etaBin;
369 // Fill the grid but do not count id in EMCal acceptance
370 //------------------------------------------------------
371 if((eta >= fEtaMin && eta <= fEtaMax) &&
372 (phi >= fPhiMin && phi <= fPhiMax)){
373 etaBin = etaBin2 = etaBin3 = -1;
379 if(eta <= fEta->At(0)) {
381 } else if(eta >= fEta->At(fNeta)) {
384 for(Int_t i=0; i<fNeta+1; i++) {
385 if(eta < fEta->At(i)) {
393 if((phi>=fPhiMin) && (phi<=fPhiMax)) {
394 if(eta <= fEta->At(0)) {
396 } else if(eta >= fEta->At(fNeta)) {
397 etaBin2 = fNeta-fEtaBinInEMCalAcc;
399 for(Int_t i=0; i<fNeta+1; i++) {
400 if(eta < fEta->At(i) && ((fEta->At(0)<eta && eta<fEtaMin) ||
401 (fEtaMax<eta && eta<fEta->At(fNeta)))) {
402 // cout << "i : " << i << endl;
405 else etaBin2 = i-1-fEtaBinInEMCalAcc;
412 if(eta <= fEta->At(0)) {
414 } else if(eta >= fEta->At(fNeta)) {
417 for(Int_t i=0; i<fNeta+1; i++) {
418 if(eta < fEta->At(i)) {
427 if(phi <= fPhi->At(0)) {
429 } else if(phi >= fPhi->At(fNphi)) {
432 for(Int_t i=0; i<fNphi+1; i++) {
433 if(phi < fPhi->At(i)) {
441 //-------------------------------------------------------
443 absID = phiBin*(fNeta+1) + etaBin;
444 if(phi>=fPhiMin && phi<=fPhiMax) {
445 if(eta >= fEtaMin && eta <= fEtaMax) absID = -1;
447 absID = (fNbinPhi+1)*(fNeta+1)
448 + (phiBin-fNbinPhi-1)*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
453 absID = (fNbinPhi+1)*(fNeta+1)
454 + fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
455 + (phiBin-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)
464 //__________________________________________________________
465 void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
467 // Get (eta,phi) position for a given index BUT loop over all entries (takes time)
469 for(Int_t j=0; j<fNphi+1; j++) {
470 for(Int_t i=0; i<fNeta+1; i++) {
473 //-------------------------------------
475 if(j*(fNeta+1)+i == index) {
482 //-------------------------------------
485 if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
486 if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
489 if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
494 if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
495 ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
496 if(ii==0) {Int_t ind = 0; eta = fEta->At(ind);}
497 else eta = fEta->At(i);
501 if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) &&
502 ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
503 +(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
512 //__________________________________________________________
513 Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta)
515 // Get index value for a (eta,phi) position - Direct value
517 Int_t ieta = GetIndexJFromEta(eta);
518 Int_t iphi = GetIndexIFromPhi(phi);
520 Int_t index = GetMatrixIndex(iphi,ieta);
523 cout << "(phi,eta) : " << phi << ", " << eta << endl;
524 cout << "index : " << index << endl;
530 //__________________________________________________________
531 Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
537 Double_t temp = (eta+fMaxEta)/(fMaxEta-fMinEta)*fNeta;
540 cout << "eta : " << eta << endl;
541 cout << "fMaxEta : " << fMaxEta << endl;
542 cout << "fMinEta : " << fMinEta << endl;
543 cout << "fNeta : " << fNeta << endl;
544 cout << "temp eta before cast : " << temp << endl;
546 idEta = static_cast<Int_t>(temp+0.5);
547 if(fDebug>20) cout << "temp eta after cast : " << idEta << endl;
550 //__________________________________________________________
551 Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
558 if(fMinPhi==0) temp = phi/(fMaxPhi-fMinPhi)*fNphi;
559 else temp = (phi-fMinPhi)/(fMaxPhi-fMinPhi)*fNphi;
563 cout << "phi : " << phi << endl;
564 cout << "fMaxPhi : " << fMaxPhi << endl;
565 cout << "fMinPhi : " << fMinPhi << endl;
566 cout << "fNphi : " << fNphi << endl;
567 cout << "temp phi before cast : " << temp << endl;
569 idPhi = static_cast<Int_t>(temp+0.5);
570 if(fDebug>20) cout << "temp phi after cast : " << idPhi << endl;
576 //__________________________________________________________
577 void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par)
579 // Allows to set parameters using only one index (if fGrid==0) !!
582 Int_t iphi = (Int_t)i/fNeta;
583 Int_t ieta = i-iphi*fNeta;
584 SetMatrixIndex(iphi,ieta,par);
589 //__________________________________________________________
590 void AliJetGrid::SetMatrixIndexes()
592 // Fill the final matrix object with the corresponding index in eta/phi
594 for(Int_t i=0; i<fNphi+1; i++){
595 for(Int_t j=0; j<fNeta+1; j++){
596 (*fIndex)(i,j) = GetIndexFromEtaPhi(fPhi->At(i),fEta->At(j))+1;
598 cout << "(*fIndex)(" << i << "," << j << ") : " << (*fIndex)(i,j) <<
599 ", phi : " << fPhi->At(i) << ", eta : " << fEta->At(j) << endl;
603 printf("##########################################################\n");
604 printf("TMatrix object filled !\n");
605 printf("Size of the object : phi x eta = (fNphi+1) x (fNeta+1) = %d\n",(fNphi+1)*(fNeta+1));
606 printf("##########################################################\n");
609 //__________________________________________________________
610 void AliJetGrid::SetIndexIJ()
612 // Set the grid index
614 for(Int_t i=0; i<fNphi+1; i++){
615 for(Int_t j=0; j<fNeta+1; j++){
616 Int_t id = static_cast<Int_t>((*fIndex)(i,j));
626 printf("##########################################################\n");
627 printf(" In SetIndexIJ - Grid indexes setted !\n");
628 printf("##########################################################\n");
631 //__________________________________________________________
632 void AliJetGrid::GetIJFromIndex(Int_t index, Int_t& i, Int_t& j) const
634 // Returns i position id of eta and j position id of phi for a given grid index
635 i = (*fIndexI)[index];
636 j = (*fIndexJ)[index];
639 //__________________________________________________________
640 void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
642 // Returns eta, phi values for a given grid index
644 phi = fPhi->At((*fIndexI)[index]);
645 eta = fEta->At((*fIndexJ)[index]);
648 //__________________________________________________________
649 Int_t AliJetGrid::GetNEntries()
651 // Returns the number of entries of the grid
654 cout << "fMaxPhi : " << fMaxPhi << endl;
655 cout << "fMaxEta : " << fMaxEta << endl;
658 Int_t indexNum = GetIndex(fMaxPhi,fMaxEta);
659 if(fDebug>20) cout << "indexNum : " << indexNum << endl;
664 //__________________________________________________________
665 Int_t AliJetGrid::GetNEntries2()
667 // Returns the number of entries of the grid
669 Int_t indexNum = GetIndex(fMaxPhi-1.,fMaxEta-0.5);
670 if(fDebug>20) cout << "indexNum : " << indexNum << endl;