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@ires.in2p3.fr
34 //=========================================================================
37 #include <Riostream.h>
44 #include "AliJetGrid.h"
49 //__________________________________________________________
50 AliJetGrid::AliJetGrid():
75 // Default constructor
78 //__________________________________________________________
79 AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax):
104 // Standard constructor
105 fPhi = new TArrayD(fNphi+1);
106 fEta = new TArrayD(fNeta+1);
107 fIndexI = new TArrayI((fNeta+1)*(fNphi+1)+1);
108 fIndexJ = new TArrayI((fNeta+1)*(fNphi+1)+1);
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;
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;
118 fIndex = new TMatrixD(fNphi+1,fNeta+1);
122 //__________________________________________________________
123 AliJetGrid::AliJetGrid(const AliJetGrid& grid) : TNamed(grid),
130 fIndexI(grid.fIndexI),
131 fIndexJ(grid.fIndexJ),
132 fPhiMin(grid.fPhiMin),
133 fPhiMax(grid.fPhiMax),
134 fEtaMin(grid.fEtaMin),
135 fEtaMax(grid.fEtaMax),
136 fEtaBinInTPCAcc(grid.fEtaBinInTPCAcc),
137 fPhiBinInTPCAcc(grid.fPhiBinInTPCAcc),
138 fEtaBinInEMCalAcc(grid.fEtaBinInEMCalAcc),
139 fPhiBinInEMCalAcc(grid.fPhiBinInEMCalAcc),
140 fNbinEta(grid.fNbinEta),
141 fNbinPhi(grid.fNbinPhi),
142 fMaxPhi(grid.fMaxPhi),
143 fMinPhi(grid.fMinPhi),
144 fMaxEta(grid.fMaxEta),
145 fMinEta(grid.fMinEta),
151 fPhi = new TArrayD(fNphi+1);
152 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = grid.fPhi->At(i);
153 fEta = new TArrayD(fNeta+1);
154 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = grid.fEta->At(i);
156 fIndex = new TMatrixD(fNphi+1,fNeta+1);
157 for(Int_t i=0; i<fNphi+1; i++) {
158 for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*grid.fIndex)(i,j);
163 AliJetGrid& AliJetGrid::operator=(const AliJetGrid& other)
172 fIndexI = other.fIndexI;
173 fIndexJ = other.fIndexJ;
174 fPhiMin = other.fPhiMin;
175 fPhiMax = other.fPhiMax;
176 fEtaMin = other.fEtaMin;
177 fEtaMax = other.fEtaMax;
178 fEtaBinInTPCAcc = other.fEtaBinInTPCAcc;
179 fPhiBinInTPCAcc = other.fPhiBinInTPCAcc;
180 fEtaBinInEMCalAcc = other.fEtaBinInEMCalAcc;
181 fPhiBinInEMCalAcc = other.fPhiBinInEMCalAcc;
182 fNbinEta = other.fNbinEta;
183 fNbinPhi = other.fNbinPhi;
184 fMaxPhi = other.fMaxPhi;
185 fMinPhi = other.fMinPhi;
186 fMaxEta = other.fMaxEta;
187 fMinEta = other.fMinEta;
188 fDebug = other.fDebug;
189 fPhi = new TArrayD(fNphi+1);
190 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = other.fPhi->At(i);
191 fEta = new TArrayD(fNeta+1);
192 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = other.fEta->At(i);
194 fIndex = new TMatrixD(fNphi+1,fNeta+1);
195 for(Int_t i=0; i<fNphi+1; i++) {
196 for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*other.fIndex)(i,j);
201 //__________________________________________________________
202 AliJetGrid::~AliJetGrid() {
212 //__________________________________________________________
213 void AliJetGrid::InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal)
215 // To set initial parameters
217 fPhiMin = phiMinCal; // rad
218 fPhiMax = phiMaxCal; // rad
221 fNbinPhi = static_cast<int>(fPhiMin/(2*TMath::Pi()/fNphi));
223 // Define some binning numbers
225 for(Int_t i=0; i<fNphi+1; i++) fPhiBinInTPCAcc++;
226 fPhiBinInEMCalAcc = 0;
228 for(Int_t i=0; i<fNeta+1; i++) fEtaBinInTPCAcc++;
229 fEtaBinInEMCalAcc = 0;
233 for(Int_t i=0; i<fNphi+1; i++) {
235 if(fPhi->At(i) >= fPhiMin &&
236 fPhi->At(i) <= fPhiMax)
239 for(Int_t i=0; i<fNeta+1; i++) {
241 if((fEta->At(i) >= fEtaMin) &&
242 (fEta->At(i) <= fEtaMax))
249 //__________________________________________________________
250 TArrayD* AliJetGrid::GetArrayEta()
252 // Returns an array with the eta points
258 //__________________________________________________________
259 TArrayD* AliJetGrid::GetArrayPhi()
261 // Returns an array with the phi points
267 //__________________________________________________________
268 TMatrixD* AliJetGrid::GetIndexObject()
270 // Returns a pointer to the matrix
276 //__________________________________________________________
277 void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi,
278 Float_t &mineta, Float_t &maxeta) const
280 // Returns EMCAL acceptance initially setted
290 //__________________________________________________________
291 void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc,
292 Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi) const
294 // Returns number of bins in TPC and EMCAL geometry
296 etabintpc = fEtaBinInTPCAcc;
297 phibintpc = fPhiBinInTPCAcc;
298 etabinemc = fEtaBinInEMCalAcc;
299 phibinemc = fPhiBinInEMCalAcc;
303 //__________________________________________________________
304 Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const
306 // Tells the index value of a corresponding (eta,phi) real position
307 // Loop over all entries -> takes time.
308 // Used one time at the begining to fill the grids
310 /* this is how bins are numbered
312 in all TPC or in TPC - EMCAL
315 ... ... . 10 | | | 11
316 ---+---+--- ---------------
317 ^ 6 | 7 | 8 or 8 | | | 9
318 | ---+---+--- ---------------
319 3 | 4 | 5 4 | 5 | 6 | 7
320 phi ---+---+--- ---------------
321 0 | 1 | 2 0 | 1 | 2 | 3
327 Int_t etaBin=0,phiBin=0,absID=0;
328 Int_t etaBin2=0,etaBin3=0;
330 // Fill all the grid in eta/phi (all TPC acceptance)
331 //-----------------------------------------------------
333 if(eta <= fEta->At(0)) {
335 } else if(eta >= fEta->At(fNeta)) {
338 for(Int_t i=0; i<fNeta+1; i++) {
339 if(eta < fEta->At(i)) {
345 if(phi <= fPhi->At(0)) {
347 } else if(phi >= fPhi->At(fNphi)) {
350 for(Int_t i=0; i<fNphi+1; i++) {
351 if(phi < fPhi->At(i)) {
358 // Calculate absolute id
359 absID = phiBin*(fNeta+1) + etaBin;
363 // Fill the grid but do not count id in EMCal acceptance
364 //------------------------------------------------------
365 if((eta >= fEtaMin && eta <= fEtaMax) &&
366 (phi >= fPhiMin && phi <= fPhiMax)){
367 etaBin = etaBin2 = etaBin3 = -1;
373 if(eta <= fEta->At(0)) {
375 } else if(eta >= fEta->At(fNeta)) {
378 for(Int_t i=0; i<fNeta+1; i++) {
379 if(eta < fEta->At(i)) {
387 if((phi>=fPhiMin) && (phi<=fPhiMax)) {
388 if(eta <= fEta->At(0)) {
390 } else if(eta >= fEta->At(fNeta)) {
391 etaBin2 = fNeta-fEtaBinInEMCalAcc;
393 for(Int_t i=0; i<fNeta+1; i++) {
394 if(eta < fEta->At(i) && ((fEta->At(0)<eta && eta<fEtaMin) ||
395 (fEtaMax<eta && eta<fEta->At(fNeta)))) {
396 // cout << "i : " << i << endl;
399 else etaBin2 = i-1-fEtaBinInEMCalAcc;
406 if(eta <= fEta->At(0)) {
408 } else if(eta >= fEta->At(fNeta)) {
411 for(Int_t i=0; i<fNeta+1; i++) {
412 if(eta < fEta->At(i)) {
421 if(phi <= fPhi->At(0)) {
423 } else if(phi >= fPhi->At(fNphi)) {
426 for(Int_t i=0; i<fNphi+1; i++) {
427 if(phi < fPhi->At(i)) {
435 //-------------------------------------------------------
437 absID = phiBin*(fNeta+1) + etaBin;
438 if(phi>=fPhiMin && phi<=fPhiMax) {
439 if(eta >= fEtaMin && eta <= fEtaMax) absID = -1;
441 absID = (fNbinPhi+1)*(fNeta+1)
442 + (phiBin-fNbinPhi-1)*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
447 absID = (fNbinPhi+1)*(fNeta+1)
448 + fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
449 + (phiBin-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)
458 //__________________________________________________________
459 void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
461 // Get (eta,phi) position for a given index BUT loop over all entries (takes time)
463 for(Int_t j=0; j<fNphi+1; j++) {
464 for(Int_t i=0; i<fNeta+1; i++) {
467 //-------------------------------------
469 if(j*(fNeta+1)+i == index) {
476 //-------------------------------------
479 if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
480 if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
483 if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
488 if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
489 ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
490 if(ii==0) {Int_t ind = 0; eta = fEta->At(ind);}
491 else eta = fEta->At(i);
495 if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) &&
496 ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
497 +(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
506 //__________________________________________________________
507 Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta)
509 // Get index value for a (eta,phi) position - Direct value
511 Int_t ieta = GetIndexJFromEta(eta);
512 Int_t iphi = GetIndexIFromPhi(phi);
514 Int_t index = GetMatrixIndex(iphi,ieta);
517 cout << "(phi,eta) : " << phi << ", " << eta << endl;
518 cout << "index : " << index << endl;
524 //__________________________________________________________
525 Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
531 Double_t temp = (eta+fMaxEta)/(fMaxEta-fMinEta)*fNeta;
534 cout << "eta : " << eta << endl;
535 cout << "fMaxEta : " << fMaxEta << endl;
536 cout << "fMinEta : " << fMinEta << endl;
537 cout << "fNeta : " << fNeta << endl;
538 cout << "temp eta before cast : " << temp << endl;
540 idEta = static_cast<Int_t>(temp+0.5);
541 if(fDebug>20) cout << "temp eta after cast : " << idEta << endl;
544 //__________________________________________________________
545 Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
552 if(fMinPhi==0) temp = phi/(fMaxPhi-fMinPhi)*fNphi;
553 else temp = (phi-fMinPhi)/(fMaxPhi-fMinPhi)*fNphi;
557 cout << "phi : " << phi << endl;
558 cout << "fMaxPhi : " << fMaxPhi << endl;
559 cout << "fMinPhi : " << fMinPhi << endl;
560 cout << "fNphi : " << fNphi << endl;
561 cout << "temp phi before cast : " << temp << endl;
563 idPhi = static_cast<Int_t>(temp+0.5);
564 if(fDebug>20) cout << "temp phi after cast : " << idPhi << endl;
570 //__________________________________________________________
571 void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par)
573 // Allows to set parameters using only one index (if fGrid==0) !!
576 Int_t iphi = (Int_t)i/fNeta;
577 Int_t ieta = i-iphi*fNeta;
578 SetMatrixIndex(iphi,ieta,par);
583 //__________________________________________________________
584 void AliJetGrid::SetMatrixIndexes()
586 // Fill the final matrix object with the corresponding index in eta/phi
588 for(Int_t i=0; i<fNphi+1; i++){
589 for(Int_t j=0; j<fNeta+1; j++){
590 (*fIndex)(i,j) = GetIndexFromEtaPhi(fPhi->At(i),fEta->At(j))+1;
592 cout << "(*fIndex)(" << i << "," << j << ") : " << (*fIndex)(i,j) <<
593 ", phi : " << fPhi->At(i) << ", eta : " << fEta->At(j) << endl;
597 printf("##########################################################\n");
598 printf("TMatrix object filled !\n");
599 printf("Size of the object : phi x eta = (fNphi+1) x (fNeta+1) = %d\n",(fNphi+1)*(fNeta+1));
600 printf("##########################################################\n");
603 //__________________________________________________________
604 void AliJetGrid::SetIndexIJ()
606 // Set the grid index
608 for(Int_t i=0; i<fNphi+1; i++){
609 for(Int_t j=0; j<fNeta+1; j++){
610 Int_t id = static_cast<Int_t>((*fIndex)(i,j));
620 printf("##########################################################\n");
621 printf(" In SetIndexIJ - Grid indexes setted !\n");
622 printf("##########################################################\n");
625 //__________________________________________________________
626 void AliJetGrid::GetIJFromIndex(Int_t index, Int_t i, Int_t j) const
628 // Returns i position id of eta and j position id of phi for a given grid index
629 i = (*fIndexI)[index];
630 j = (*fIndexJ)[index];
633 //__________________________________________________________
634 void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
636 // Returns eta, phi values for a given grid index
638 phi = fPhi->At((*fIndexI)[index]);
639 eta = fEta->At((*fIndexJ)[index]);
642 //__________________________________________________________
643 Int_t AliJetGrid::GetNEntries()
645 // Returns the number of entries of the grid
648 cout << "fMaxPhi : " << fMaxPhi << endl;
649 cout << "fMaxEta : " << fMaxEta << endl;
652 Int_t indexNum = GetIndex(fMaxPhi,fMaxEta);
653 if(fDebug>20) cout << "indexNum : " << indexNum << endl;
658 //__________________________________________________________
659 Int_t AliJetGrid::GetNEntries2()
661 // Returns the number of entries of the grid
663 Int_t indexNum = GetIndex(fMaxPhi-1.,fMaxEta-0.5);
664 if(fDebug>20) cout << "indexNum : " << indexNum << endl;