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"
50 //__________________________________________________________
51 AliJetGrid::AliJetGrid():
76 // Default constructor
79 //__________________________________________________________
80 AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax):
105 // Standard constructor
106 fPhi = new TArrayD(fNphi+1);
107 fEta = new TArrayD(fNeta+1);
108 fIndexI = new TArrayI((fNeta+1)*(fNphi+1)+1);
109 fIndexJ = new TArrayI((fNeta+1)*(fNphi+1)+1);
111 for(Int_t i=0; i<fNphi+1; i++) {
112 if(fNphi!=0) (*fPhi)[i] = (phiMax-phiMin)/fNphi*i+phiMin;
113 else (*fPhi)[i] = phiMin+(phiMax-phiMin)/2;
115 for(Int_t i=0; i<fNeta+1; i++) {
116 if(fNeta!=0) (*fEta)[i] = (etaMax-etaMin)/fNeta*i+etaMin;
117 else (*fEta)[i] = etaMin+(etaMax-etaMin)/2;
121 for(Int_t i=0; i<(fNphi+1); i++) cout << (*fPhi)[i] << endl;
122 for(Int_t i=0; i<(fNeta+1); i++) cout << (*fEta)[i] << endl;
125 fIndex = new TMatrixD(fNphi+1,fNeta+1);
129 //__________________________________________________________
130 AliJetGrid::AliJetGrid(const AliJetGrid& grid) :
138 fIndexI(grid.fIndexI),
139 fIndexJ(grid.fIndexJ),
140 fPhiMin(grid.fPhiMin),
141 fPhiMax(grid.fPhiMax),
142 fEtaMin(grid.fEtaMin),
143 fEtaMax(grid.fEtaMax),
144 fEtaBinInTPCAcc(grid.fEtaBinInTPCAcc),
145 fPhiBinInTPCAcc(grid.fPhiBinInTPCAcc),
146 fEtaBinInEMCalAcc(grid.fEtaBinInEMCalAcc),
147 fPhiBinInEMCalAcc(grid.fPhiBinInEMCalAcc),
148 fNbinEta(grid.fNbinEta),
149 fNbinPhi(grid.fNbinPhi),
150 fMaxPhi(grid.fMaxPhi),
151 fMinPhi(grid.fMinPhi),
152 fMaxEta(grid.fMaxEta),
153 fMinEta(grid.fMinEta),
159 fPhi = new TArrayD(fNphi+1);
160 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = grid.fPhi->At(i);
161 fEta = new TArrayD(fNeta+1);
162 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = grid.fEta->At(i);
164 fIndex = new TMatrixD(fNphi+1,fNeta+1);
165 for(Int_t i=0; i<fNphi+1; i++) {
166 for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*grid.fIndex)(i,j);
171 AliJetGrid& AliJetGrid::operator=(const AliJetGrid& other)
174 if (this != &other) {
181 fIndexI = other.fIndexI;
182 fIndexJ = other.fIndexJ;
183 fPhiMin = other.fPhiMin;
184 fPhiMax = other.fPhiMax;
185 fEtaMin = other.fEtaMin;
186 fEtaMax = other.fEtaMax;
187 fEtaBinInTPCAcc = other.fEtaBinInTPCAcc;
188 fPhiBinInTPCAcc = other.fPhiBinInTPCAcc;
189 fEtaBinInEMCalAcc = other.fEtaBinInEMCalAcc;
190 fPhiBinInEMCalAcc = other.fPhiBinInEMCalAcc;
191 fNbinEta = other.fNbinEta;
192 fNbinPhi = other.fNbinPhi;
193 fMaxPhi = other.fMaxPhi;
194 fMinPhi = other.fMinPhi;
195 fMaxEta = other.fMaxEta;
196 fMinEta = other.fMinEta;
197 fDebug = other.fDebug;
198 fPhi = new TArrayD(fNphi+1);
199 for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = other.fPhi->At(i);
200 fEta = new TArrayD(fNeta+1);
201 for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = other.fEta->At(i);
203 fIndex = new TMatrixD(fNphi+1,fNeta+1);
204 for(Int_t i=0; i<fNphi+1; i++) {
205 for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*other.fIndex)(i,j);
212 //__________________________________________________________
213 AliJetGrid::~AliJetGrid() {
223 //__________________________________________________________
224 void AliJetGrid::InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal)
226 // To set initial parameters
228 fPhiMin = phiMinCal; // rad
229 fPhiMax = phiMaxCal; // rad
232 fNbinPhi = static_cast<int>(fPhiMin/(2*TMath::Pi()/fNphi));
234 // Define some binning numbers
236 for(Int_t i=0; i<fNphi+1; i++) fPhiBinInTPCAcc++;
237 fPhiBinInEMCalAcc = 0;
239 for(Int_t i=0; i<fNeta+1; i++) fEtaBinInTPCAcc++;
240 fEtaBinInEMCalAcc = 0;
244 for(Int_t i=0; i<fNphi+1; i++) {
246 if(fPhi->At(i) >= fPhiMin &&
247 fPhi->At(i) <= fPhiMax)
250 for(Int_t i=0; i<fNeta+1; i++) {
252 if((fEta->At(i) >= fEtaMin) &&
253 (fEta->At(i) <= fEtaMax))
260 //__________________________________________________________
261 TArrayD* AliJetGrid::GetArrayEta()
263 // Returns an array with the eta points
269 //__________________________________________________________
270 TArrayD* AliJetGrid::GetArrayPhi()
272 // Returns an array with the phi points
278 //__________________________________________________________
279 TMatrixD* AliJetGrid::GetIndexObject()
281 // Returns a pointer to the matrix
287 //__________________________________________________________
288 void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi,
289 Float_t &mineta, Float_t &maxeta) const
291 // Returns EMCAL acceptance initially setted
301 //__________________________________________________________
302 void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc,
303 Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi) const
305 // Returns number of bins in TPC and EMCAL geometry
307 etabintpc = fEtaBinInTPCAcc;
308 phibintpc = fPhiBinInTPCAcc;
309 etabinemc = fEtaBinInEMCalAcc;
310 phibinemc = fPhiBinInEMCalAcc;
314 //__________________________________________________________
315 Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const
317 // Tells the index value of a corresponding (eta,phi) real position
318 // Loop over all entries -> takes time.
319 // Used one time at the begining to fill the grids
321 /* this is how bins are numbered
323 in all TPC or in TPC - EMCAL
326 ... ... . 10 | | | 11
327 ---+---+--- ---------------
328 ^ 6 | 7 | 8 or 8 | | | 9
329 | ---+---+--- ---------------
330 3 | 4 | 5 4 | 5 | 6 | 7
331 phi ---+---+--- ---------------
332 0 | 1 | 2 0 | 1 | 2 | 3
338 Int_t etaBin=0,phiBin=0,absID=0;
339 Int_t etaBin2=0,etaBin3=0;
341 // Fill all the grid in eta/phi (all TPC acceptance)
342 //-----------------------------------------------------
344 if(eta <= fEta->At(0)) {
346 } else if(eta >= fEta->At(fNeta)) {
349 for(Int_t i=0; i<fNeta+1; i++) {
350 if(eta < fEta->At(i)) {
356 if(phi <= fPhi->At(0)) {
358 } else if(phi >= fPhi->At(fNphi)) {
361 for(Int_t i=0; i<fNphi+1; i++) {
362 if(phi < fPhi->At(i)) {
369 // Calculate absolute id
370 absID = phiBin*(fNeta+1) + etaBin;
374 // Fill the grid but do not count id in EMCal acceptance
375 //------------------------------------------------------
376 if((eta >= fEtaMin && eta <= fEtaMax) &&
377 (phi >= fPhiMin && phi <= fPhiMax)){
378 etaBin = etaBin2 = etaBin3 = -1;
384 if(eta <= fEta->At(0)) {
386 } else if(eta >= fEta->At(fNeta)) {
389 for(Int_t i=0; i<fNeta+1; i++) {
390 if(eta < fEta->At(i)) {
398 if((phi>=fPhiMin) && (phi<=fPhiMax)) {
399 if(eta <= fEta->At(0)) {
401 } else if(eta >= fEta->At(fNeta)) {
402 etaBin2 = fNeta-fEtaBinInEMCalAcc;
404 for(Int_t i=0; i<fNeta+1; i++) {
405 if(eta < fEta->At(i) && ((fEta->At(0)<eta && eta<fEtaMin) ||
406 (fEtaMax<eta && eta<fEta->At(fNeta)))) {
407 // cout << "i : " << i << endl;
410 else etaBin2 = i-1-fEtaBinInEMCalAcc;
417 if(eta <= fEta->At(0)) {
419 } else if(eta >= fEta->At(fNeta)) {
422 for(Int_t i=0; i<fNeta+1; i++) {
423 if(eta < fEta->At(i)) {
432 if(phi <= fPhi->At(0)) {
434 } else if(phi >= fPhi->At(fNphi)) {
437 for(Int_t i=0; i<fNphi+1; i++) {
438 if(phi < fPhi->At(i)) {
446 //-------------------------------------------------------
448 absID = phiBin*(fNeta+1) + etaBin;
449 if(phi>=fPhiMin && phi<=fPhiMax) {
450 if(eta >= fEtaMin && eta <= fEtaMax) absID = -1;
452 absID = (fNbinPhi+1)*(fNeta+1)
453 + (phiBin-fNbinPhi-1)*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
458 absID = (fNbinPhi+1)*(fNeta+1)
459 + fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
460 + (phiBin-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)
469 //__________________________________________________________
470 void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
472 // Get (eta,phi) position for a given index BUT loop over all entries (takes time)
474 for(Int_t j=0; j<fNphi+1; j++) {
475 for(Int_t i=0; i<fNeta+1; i++) {
478 //-------------------------------------
480 if(j*(fNeta+1)+i == index) {
487 //-------------------------------------
490 if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
491 if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
494 if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
499 if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
500 ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
501 if(ii==0) {Int_t ind = 0; eta = fEta->At(ind);}
502 else eta = fEta->At(i);
506 if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) &&
507 ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
508 +(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
517 //__________________________________________________________
518 Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta)
520 // Get index value for a (eta,phi) position - Direct value
522 Int_t ieta = GetIndexJFromEta(eta);
523 Int_t iphi = GetIndexIFromPhi(phi);
525 Int_t index = GetMatrixIndex(iphi,ieta);
528 cout << "(phi,eta) : " << phi << ", " << eta << endl;
529 cout << "index : " << index << endl;
535 //__________________________________________________________
536 Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
542 Double_t temp = (eta+fMaxEta)/(fMaxEta-fMinEta)*fNeta;
545 cout << "eta : " << eta << endl;
546 cout << "fMaxEta : " << fMaxEta << endl;
547 cout << "fMinEta : " << fMinEta << endl;
548 cout << "fNeta : " << fNeta << endl;
549 cout << "temp eta before cast : " << temp << endl;
551 idEta = static_cast<Int_t>(temp+0.5);
552 if(fDebug>20) cout << "temp eta after cast : " << idEta << endl;
555 //__________________________________________________________
556 Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
563 if(fMinPhi==0) temp = phi/(fMaxPhi-fMinPhi)*fNphi;
564 else temp = (phi-fMinPhi)/(fMaxPhi-fMinPhi)*fNphi;
568 cout << "phi : " << phi << endl;
569 cout << "fMaxPhi : " << fMaxPhi << endl;
570 cout << "fMinPhi : " << fMinPhi << endl;
571 cout << "fNphi : " << fNphi << endl;
572 cout << "temp phi before cast : " << temp << endl;
574 idPhi = static_cast<Int_t>(temp+0.5);
575 if(fDebug>20) cout << "temp phi after cast : " << idPhi << endl;
581 //__________________________________________________________
582 void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par)
584 // Allows to set parameters using only one index (if fGrid==0) !!
587 Int_t iphi = (Int_t)i/fNeta;
588 Int_t ieta = i-iphi*fNeta;
589 SetMatrixIndex(iphi,ieta,par);
594 //__________________________________________________________
595 void AliJetGrid::SetMatrixIndexes()
597 // Fill the final matrix object with the corresponding index in eta/phi
599 for(Int_t i=0; i<fNphi+1; i++){
600 for(Int_t j=0; j<fNeta+1; j++){
601 (*fIndex)(i,j) = GetIndexFromEtaPhi(fPhi->At(i),fEta->At(j))+1;
603 cout << "(*fIndex)(" << i << "," << j << ") : " << (*fIndex)(i,j) <<
604 ", phi : " << fPhi->At(i) << ", eta : " << fEta->At(j) << endl;
608 printf("##########################################################\n");
609 printf("TMatrix object filled !\n");
610 printf("Size of the object : phi x eta = (fNphi+1) x (fNeta+1) = %d\n",(fNphi+1)*(fNeta+1));
611 printf("##########################################################\n");
614 //__________________________________________________________
615 void AliJetGrid::SetIndexIJ()
617 // Set the grid index
619 for(Int_t i=0; i<fNphi+1; i++){
620 for(Int_t j=0; j<fNeta+1; j++){
621 Int_t id = static_cast<Int_t>((*fIndex)(i,j));
631 printf("##########################################################\n");
632 printf(" In SetIndexIJ - Grid indexes setted !\n");
633 printf("##########################################################\n");
636 //__________________________________________________________
637 void AliJetGrid::GetIJFromIndex(Int_t index, Int_t& i, Int_t& j) const
639 // Returns i position id of eta and j position id of phi for a given grid index
640 i = (*fIndexI)[index];
641 j = (*fIndexJ)[index];
644 //__________________________________________________________
645 void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
647 // Returns eta, phi values for a given grid index
649 phi = fPhi->At((*fIndexI)[index]);
650 eta = fEta->At((*fIndexJ)[index]);
653 //__________________________________________________________
654 Int_t AliJetGrid::GetNEntries()
656 // Returns the number of entries of the grid
659 cout << "fMaxPhi : " << fMaxPhi << endl;
660 cout << "fMaxEta : " << fMaxEta << endl;
663 Int_t indexNum = GetIndex(fMaxPhi,fMaxEta);
664 if(fDebug>20) cout << "indexNum : " << indexNum << endl;
669 //__________________________________________________________
670 Int_t AliJetGrid::GetNEntries2()
672 // Returns the number of entries of the grid
674 Int_t indexNum = GetIndex(fMaxPhi-1.,fMaxEta-0.5);
675 if(fDebug>20) cout << "indexNum : " << indexNum << endl;