JETAN/AliJetGrid.cxx

   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
  39 #include <TMatrixD.h>
  40 #include <TArrayD.h>
  41 #include <TArrayI.h>
  42 // AliRoot headers
  43 #include "AliJetGrid.h"
  44
  45 ClassImp(AliJetGrid)
  46
  47 //__________________________________________________________
  48 AliJetGrid::AliJetGrid() {
  49
  50   // Default constructor
  51
  52   fNphi = 0;
  53   fNeta = 0;
  54   fPhi = 0;
  55   fEta = 0;
  56   fIndex = 0x0;
  57   fPhiMin = 0.; // acceptance removed inside
  58   fPhiMax = 0.; // acceptance removed inside
  59   fEtaMin = 0.; // acceptance removed inside
  60   fEtaMax = 0.; // acceptance removed inside
  61   fMinPhi = 0.; // total acceptance
  62   fMaxPhi = 0.; // total acceptance
  63   fMinEta = 0.; // total acceptance
  64   fMaxEta = 0.; // total acceptance
  65   fEtaBinInTPCAcc = 0;
  66   fPhiBinInTPCAcc = 0;
  67   fEtaBinInEMCalAcc = 0;
  68   fPhiBinInEMCalAcc = 0;
  69   fNbinPhi = 0;
  70
  71   fDebug = 1;
  72 }
  73
  74 //__________________________________________________________
  75 AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax) {
  76
  77   // Standard constructor
  78
  79   fNphi = nphi;
  80   fNeta = neta;
  81   fPhiMin = 0.; // rad - acceptance removed inside
  82   fPhiMax = 0.; // rad - acceptance removed inside
  83   fEtaMin = 0.; //  acceptance removed inside
  84   fEtaMax = 0.; //  acceptance removed inside
  85   fNbinPhi = 0;
  86
  87   fMaxPhi = phiMax; // total acceptance
  88   fMinPhi = phiMin; // total acceptance
  89   fMaxEta = etaMax; // total acceptance
  90   fMinEta = etaMin; // total acceptance
  91
  92   fPhi  = new TArrayD(fNphi+1);
  93   fEta  = new TArrayD(fNeta+1);
  94   fIndexI = new TArrayI((fNeta+1)*(fNphi+1)+1);
  95   fIndexJ = new TArrayI((fNeta+1)*(fNphi+1)+1);
  96
  97   for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = (phiMax-phiMin)/fNphi*i+phiMin;
  98   for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = (etaMax-etaMin)/fNeta*i+etaMin;
  99
 100   if(fDebug > 3){
 101     for(Int_t i=0; i<(fNphi+1); i++)  cout << (*fPhi)[i] << endl;
 102     for(Int_t i=0; i<(fNeta+1); i++)  cout << (*fEta)[i] << endl;
 103   }
 104
 105   fIndex = new TMatrixD(fNphi+1,fNeta+1);
 106
 107 }
 108
 109 //__________________________________________________________
 110 AliJetGrid::AliJetGrid(const AliJetGrid& grid):TNamed(grid) {
 111
 112   // Copy constructor
 113
 114   fNphi = grid.fNphi;
 115   fNeta = grid.fNeta;
 116   fPhiMin = grid.fPhiMin;
 117   fPhiMax = grid.fPhiMax;
 118   fEtaMin = grid.fEtaMin;
 119   fEtaMax = grid.fEtaMax;
 120   fEtaBinInTPCAcc = grid.fEtaBinInTPCAcc;
 121   fPhiBinInTPCAcc = grid.fPhiBinInTPCAcc;
 122   fEtaBinInEMCalAcc = grid.fEtaBinInEMCalAcc;
 123   fPhiBinInEMCalAcc = grid.fPhiBinInEMCalAcc;
 124   fNbinPhi = grid.fNbinPhi;
 125   fMaxPhi = grid.fMaxPhi;
 126   fMinPhi = grid.fMinPhi;
 127   fMaxEta = grid.fMaxEta;
 128   fMinEta = grid.fMinEta;
 129
 130   fPhi = new TArrayD(fNphi+1);
 131   for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = grid.fPhi->At(i);
 132   fEta = new TArrayD(fNeta+1);
 133   for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = grid.fEta->At(i);
 134
 135   fIndex = new TMatrixD(fNphi+1,fNeta+1);
 136   for(Int_t i=0; i<fNphi+1; i++) {
 137     for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*grid.fIndex)(i,j);
 138   }
 139 }
 140
 141 //__________________________________________________________
 142 AliJetGrid::~AliJetGrid() {
 143
 144   // Destructor
 145   delete fPhi;
 146   delete fEta;
 147   delete fIndexI;
 148   delete fIndexJ;
 149   delete fIndex;
 150 }
 151
 152 //__________________________________________________________
 153 void AliJetGrid::InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal)
 154 { // To set initial parameters
 155
 156   fPhiMin = phiMinCal; // rad
 157   fPhiMax = phiMaxCal; // rad
 158   fEtaMin = etaMinCal;
 159   fEtaMax = etaMaxCal;
 160   fNbinPhi = static_cast<int>(fPhiMin/(2*TMath::Pi()/fNphi));
 161
 162   // Define some binning numbers
 163   if(fGrid==0){
 164     for(Int_t i=0; i<fNphi+1; i++) fPhiBinInTPCAcc++;
 165     fPhiBinInEMCalAcc = 0;
 166
 167     for(Int_t i=0; i<fNeta+1; i++) fEtaBinInTPCAcc++;
 168     fEtaBinInEMCalAcc = 0;
 169   }
 170
 171   if(fGrid==1){
 172     for(Int_t i=0; i<fNphi+1; i++) {
 173       fPhiBinInTPCAcc++;
 174       if(fPhi->At(i) >= fPhiMin &&
 175          fPhi->At(i) <= fPhiMax)
 176         fPhiBinInEMCalAcc++;
 177     }
 178     for(Int_t i=0; i<fNeta+1; i++) {
 179       fEtaBinInTPCAcc++;
 180       if((fEta->At(i) >= fEtaMin) &&
 181          (fEta->At(i) <= fEtaMax))
 182         fEtaBinInEMCalAcc++;
 183     }
 184   }
 185
 186 }
 187
 188 //__________________________________________________________
 189 TArrayD* AliJetGrid::GetArrayEta()
 190 { // Returns an array with the eta points
 191
 192   return fEta;
 193
 194 }
 195
 196 //__________________________________________________________
 197 TArrayD* AliJetGrid::GetArrayPhi()
 198 { // Returns an array with the phi points
 199
 200   return fPhi;
 201
 202 }
 203
 204 //__________________________________________________________
 205 TMatrixD* AliJetGrid::GetIndexObject()
 206 { // Returns a pointer to the matrix
 207
 208   return fIndex;
 209
 210 }
 211
 212 //__________________________________________________________
 213 void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi,
 214                                 Float_t &mineta, Float_t &maxeta)
 215 { // Returns EMCAL acceptance initially setted
 216
 217   nphi = fNphi;
 218   neta = fNeta;
 219   minphi = fPhiMin;
 220   maxphi = fPhiMax;
 221   mineta = fEtaMin;
 222   maxeta = fEtaMax;
 223 }
 224
 225 //__________________________________________________________
 226 void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc,
 227                                 Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi)
 228 { // Returns number of bins in TPC and EMCAL geometry
 229
 230   etabintpc = fEtaBinInTPCAcc;
 231   phibintpc = fPhiBinInTPCAcc;
 232   etabinemc = fEtaBinInEMCalAcc;
 233   phibinemc = fPhiBinInEMCalAcc;
 234   nbinphi = fNbinPhi;
 235 }
 236
 237 //__________________________________________________________
 238 Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const
 239 { // Tells the index value of a corresponding (eta,phi) real position
 240   // Loop over all entries -> takes time.
 241   // Used one time at the begining to fill the grids
 242
 243   /*   this is how bins are numbered
 244
 245        in all TPC        or      in TPC - EMCAL
 246                                 ... ... ... ..
 247                                 ---------------
 248           ...  ... .            10 |   |   | 11
 249           ---+---+---           ---------------
 250     ^      6 | 7 | 8     or      8 |   |   | 9
 251     |     ---+---+---           ---------------
 252            3 | 4 | 5             4 | 5 | 6 | 7
 253    phi    ---+---+---           ---------------
 254            0 | 1 | 2             0 | 1 | 2 | 3
 255
 256
 257              eta ->
 258   */
 259
 260   Int_t etaBin=0,phiBin=0,absID=0;
 261   Int_t etaBin2=0,etaBin3=0;
 262
 263   // Fill all the grid in eta/phi (all TPC acceptance)
 264   //-----------------------------------------------------
 265   if(fGrid==0){
 266     if(eta <= fEta->At(0)) {
 267       etaBin = 0;
 268     } else if(eta >= fEta->At(fNeta)) {
 269       etaBin = fNeta;
 270     } else {
 271       for(Int_t i=0; i<fNeta+1; i++) {
 272         if(eta < fEta->At(i)) {
 273           etaBin = i-1;
 274           break;
 275         }
 276       }
 277     }
 278     if(phi <= fPhi->At(0)) {
 279       phiBin = 0;
 280     } else if(phi >= fPhi->At(fNphi)) {
 281       phiBin = fNphi;
 282     } else {
 283       for(Int_t i=0; i<fNphi+1; i++) {
 284         if(phi < fPhi->At(i)) {
 285           phiBin = i-1;
 286           break;
 287         }
 288       }
 289     }
 290
 291     // Calculate absolute id
 292     absID = phiBin*(fNeta+1) + etaBin;
 293
 294   }
 295
 296   // Fill the grid but do not count id in EMCal acceptance
 297   //------------------------------------------------------
 298   if((eta >= fEtaMin && eta <= fEtaMax) &&
 299      (phi >= fPhiMin && phi <= fPhiMax)){
 300     etaBin = etaBin2 = etaBin3 = -1;
 301     phiBin = -1;
 302   }
 303
 304   if(fGrid == 1){
 305     if(phi<fPhiMin) {
 306       if(eta <= fEta->At(0)) {
 307         etaBin = 0;
 308       } else if(eta >= fEta->At(fNeta)) {
 309         etaBin = fNeta;
 310       } else {
 311         for(Int_t i=0; i<fNeta+1; i++) {
 312           if(eta < fEta->At(i)) {
 313             etaBin = i-1;
 314             break;
 315           }
 316         }
 317       }
 318     }
 319     else {
 320       if((phi>=fPhiMin) && (phi<=fPhiMax)) {
 321         if(eta <= fEta->At(0)) {
 322           etaBin2 = 0;
 323         } else if(eta >= fEta->At(fNeta)) {
 324           etaBin2 = fNeta-fEtaBinInEMCalAcc;
 325         } else {
 326           for(Int_t i=0; i<fNeta+1; i++) {
 327             if(eta < fEta->At(i) && ((fEta->At(0)<eta && eta<fEtaMin) ||
 328                                      (fEtaMax<eta && eta<fEta->At(fNeta)))) {
 329               //                 cout << "i : " << i << endl;
 330               if(eta<fEtaMin)
 331                 etaBin2 = i-1;
 332               else etaBin2 = i-1-fEtaBinInEMCalAcc;
 333               break;
 334             }
 335           }
 336         }
 337       }
 338       else {
 339         if(eta <= fEta->At(0)) {
 340           etaBin3 = 0;
 341         } else if(eta >= fEta->At(fNeta)) {
 342           etaBin3 = fNeta;
 343         } else {
 344           for(Int_t i=0; i<fNeta+1; i++) {
 345             if(eta < fEta->At(i)) {
 346               etaBin3 = i-1;
 347               break;
 348             }
 349           }
 350         }
 351       }
 352     }
 353
 354     if(phi <= fPhi->At(0)) {
 355       phiBin = 0;
 356     } else if(phi >= fPhi->At(fNphi)) {
 357       phiBin = fNphi;
 358     } else {
 359       for(Int_t i=0; i<fNphi+1; i++) {
 360         if(phi < fPhi->At(i)) {
 361           phiBin = i-1;
 362           break;
 363         }
 364       }
 365     }
 366
 367     // Calculate absID
 368     //-------------------------------------------------------
 369     if(phi<fPhiMin)
 370       absID = phiBin*(fNeta+1) + etaBin;
 371     if(phi>=fPhiMin && phi<=fPhiMax) {
 372       if(eta >= fEtaMin && eta <= fEtaMax) absID = -1;
 373       else{
 374         absID = (fNbinPhi+1)*(fNeta+1)
 375           + (phiBin-fNbinPhi-1)*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
 376           + etaBin2;
 377       }
 378     }
 379     if(phi>fPhiMax)
 380       absID = (fNbinPhi+1)*(fNeta+1)
 381         + fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
 382         + (phiBin-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)
 383         + etaBin3;
 384
 385   } // END OPTION==1
 386
 387   return absID;
 388
 389 }
 390
 391 //__________________________________________________________
 392 void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
 393 { // Get (eta,phi) position for a given index BUT loop over all entries (takes time)
 394
 395   for(Int_t j=0; j<fNphi+1; j++) {
 396     for(Int_t i=0; i<fNeta+1; i++) {
 397
 398       // TPC grid only
 399       //-------------------------------------
 400       if(fGrid==0) {
 401         if(j*(fNeta+1)+i == index) {
 402           eta = fEta->At(i);
 403           phi = fPhi->At(j);
 404         }
 405       }
 406
 407       // TPC-EMCAL grid
 408       //-------------------------------------
 409       Int_t ii = 0;
 410       if(i==0) ii = 0;
 411       if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
 412       if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
 413
 414       if(fGrid==1) {
 415         if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
 416           eta = fEta->At(i);
 417           phi = fPhi->At(j);
 418         }
 419
 420         if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
 421            ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
 422           if(ii==0) {Int_t ind = 0; eta = fEta->At(ind);}
 423           else eta = fEta->At(i);
 424           phi = fPhi->At(j);
 425         }
 426
 427         if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) &&
 428            ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
 429             +(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
 430           eta = fEta->At(i);
 431           phi = fPhi->At(j);
 432         }
 433       }
 434     }
 435   }
 436 }
 437
 438 //__________________________________________________________
 439 Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta)
 440 { // Get index value for a (eta,phi) position - Direct value
 441
 442   Int_t ieta = GetIndexJFromEta(eta);
 443   Int_t iphi = GetIndexIFromPhi(phi);
 444
 445   Int_t index = GetMatrixIndex(iphi,ieta);
 446
 447   if(fDebug>10){
 448     cout << "(phi,eta) : " << phi << ", " << eta << endl;
 449     cout << "index : " << index << endl;
 450   }
 451   return index;
 452
 453 }
 454
 455 //__________________________________________________________
 456 Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
 457 { // Get eta id
 458   // Eta discretized
 459
 460   Int_t idEta =0;
 461   Double_t temp = (eta+fMaxEta)/(fMaxEta-fMinEta)*fNeta;
 462   if(fDebug>20)
 463     {
 464       cout << "eta : " << eta << endl;
 465       cout << "fMaxEta : " << fMaxEta << endl;
 466       cout << "fMinEta : " << fMinEta << endl;
 467       cout << "fNeta : " << fNeta << endl;
 468       cout << "temp eta before cast : " << temp << endl;
 469     }
 470   idEta = static_cast<Int_t>(temp+0.5);
 471   if(fDebug>20) cout << "temp eta after cast : " << idEta << endl;
 472   return idEta;
 473 }
 474 //__________________________________________________________
 475 Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
 476 { // Get phi id
 477   // Phi discretized
 478
 479   Int_t idPhi = 0;
 480   Double_t temp = 0.;
 481   if(fMinPhi==0) temp = phi/(fMaxPhi-fMinPhi)*fNphi;
 482   else temp = (phi-fMinPhi)/(fMaxPhi-fMinPhi)*fNphi;
 483
 484   if(fDebug>20)
 485     {
 486       cout << "phi : " << phi << endl;
 487       cout << "fMaxPhi : " << fMaxPhi << endl;
 488       cout << "fMinPhi : " << fMinPhi << endl;
 489       cout << "fNphi : " << fNphi << endl;
 490       cout << "temp phi before cast : " << temp << endl;
 491     }
 492   idPhi = static_cast<Int_t>(temp+0.5);
 493   if(fDebug>20) cout << "temp phi after cast : " << idPhi << endl;
 494   return idPhi;
 495
 496
 497 }
 498
 499 //__________________________________________________________
 500 void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par)
 501 { // Allows to set parameters using only one index (if fGrid==0) !!
 502   // Not used !
 503
 504   Int_t iphi = (Int_t)i/fNeta;
 505   Int_t ieta = i-iphi*fNeta;
 506   SetMatrixIndex(iphi,ieta,par);
 507
 508   return;
 509 }
 510
 511 //__________________________________________________________
 512 void AliJetGrid::SetMatrixIndexes()
 513 { // Fill the final matrix object with the corresponding index in eta/phi
 514
 515   for(Int_t i=0; i<fNphi+1; i++){
 516     for(Int_t j=0; j<fNeta+1; j++){
 517       (*fIndex)(i,j) = GetIndexFromEtaPhi(fPhi->At(i),fEta->At(j))+1;
 518       if(fDebug>2){
 519         cout << "(*fIndex)(" << i << "," << j << ") : " << (*fIndex)(i,j) <<
 520           ", phi : " << fPhi->At(i) << ", eta : " << fEta->At(j) << endl;
 521       }
 522     }
 523   }
 524   printf("##########################################################\n");
 525   printf("TMatrix object filled !\n");
 526   printf("Size of the object : phi x eta = (fNphi+1) x (fNeta+1) = %d\n",(fNphi+1)*(fNeta+1));
 527   printf("##########################################################\n");
 528 }
 529
 530 //__________________________________________________________
 531 void AliJetGrid::SetIndexIJ()
 532 { //
 533
 534   for(Int_t i=0; i<fNphi+1; i++){
 535     for(Int_t j=0; j<fNeta+1; j++){
 536       Int_t id = static_cast<Int_t>((*fIndex)(i,j));
 537
 538       if(id!=-1)
 539         {
 540           (*fIndexI)[id] = i;
 541           (*fIndexJ)[id] = j;
 542         }
 543     }
 544   }
 545
 546   printf("##########################################################\n");
 547   printf("     In SetIndexIJ - Grid indexes setted !\n");
 548   printf("##########################################################\n");
 549 }
 550
 551 //__________________________________________________________
 552 void AliJetGrid::GetIJFromIndex(Int_t index, Int_t i, Int_t j)
 553 { // returns i position id of eta and j position id of phi for a given grid index
 554   i = (*fIndexI)[index];
 555   j = (*fIndexJ)[index];
 556 }
 557
 558 //__________________________________________________________
 559 void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
 560 { // returns eta, phi values for a given grid index
 561
 562   phi = fPhi->At((*fIndexI)[index]);
 563   eta = fEta->At((*fIndexJ)[index]);
 564 }
 565
 566 //__________________________________________________________
 567 Int_t AliJetGrid::GetNEntries()
 568 { // Returns the number of entries of the grid
 569
 570   if(fDebug>20){
 571     cout << "fMaxPhi : " << fMaxPhi << endl;
 572     cout << "fMaxEta : " << fMaxEta << endl;
 573   }
 574
 575   Int_t indexNum = GetIndex(fMaxPhi,fMaxEta);
 576   if(fDebug>20) cout << "indexNum : " << indexNum << endl;
 577   return indexNum;
 578
 579 }
 580
 581 //__________________________________________________________
 582 Int_t AliJetGrid::GetNEntries2()
 583 { // Returns the number of entries of the grid
 584
 585   Int_t indexNum = GetIndex(fMaxPhi-1.,fMaxEta-0.5);
 586     if(fDebug>20) cout << "indexNum : " << indexNum << endl;
 587   return indexNum;
 588
 589 }
 590
 591
 592
 593
 594
 595