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 <TMath.h>
  40 #include <TMatrixD.h>
  41 #include <TArrayD.h>
  42 #include <TArrayI.h>
  43 // AliRoot headers
  44 #include "AliJetGrid.h"
  45
  46 ClassImp(AliJetGrid)
  47
  48
  49 //__________________________________________________________
  50 AliJetGrid::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  {
  75   // Default constructor
  76 }
  77
  78 //__________________________________________________________
  79 AliJetGrid::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  {
 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);
 109
 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;
 112
 113   if(fDebug > 3){
 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;
 116   }
 117
 118   fIndex = new TMatrixD(fNphi+1,fNeta+1);
 119
 120 }
 121
 122 //__________________________________________________________
 123 AliJetGrid::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 //__________________________________________________________
 155 AliJetGrid::~AliJetGrid() {
 156
 157   // Destructor
 158   delete fPhi;
 159   delete fEta;
 160   delete fIndexI;
 161   delete fIndexJ;
 162   delete fIndex;
 163 }
 164
 165 //__________________________________________________________
 166 void AliJetGrid::InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal)
 167 {
 168 // To set initial parameters
 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 //__________________________________________________________
 203 TArrayD* AliJetGrid::GetArrayEta()
 204 {
 205 // Returns an array with the eta points
 206
 207   return fEta;
 208
 209 }
 210
 211 //__________________________________________________________
 212 TArrayD* AliJetGrid::GetArrayPhi()
 213 {
 214 // Returns an array with the phi points
 215
 216   return fPhi;
 217
 218 }
 219
 220 //__________________________________________________________
 221 TMatrixD* AliJetGrid::GetIndexObject()
 222 {
 223 // Returns a pointer to the matrix
 224
 225   return fIndex;
 226
 227 }
 228
 229 //__________________________________________________________
 230 void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi,
 231                                 Float_t &mineta, Float_t &maxeta) const
 232 {
 233 // Returns EMCAL acceptance initially setted
 234
 235   nphi = fNphi;
 236   neta = fNeta;
 237   minphi = fPhiMin;
 238   maxphi = fPhiMax;
 239   mineta = fEtaMin;
 240   maxeta = fEtaMax;
 241 }
 242
 243 //__________________________________________________________
 244 void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc,
 245                                 Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi) const
 246 {
 247 // Returns number of bins in TPC and EMCAL geometry
 248
 249   etabintpc = fEtaBinInTPCAcc;
 250   phibintpc = fPhiBinInTPCAcc;
 251   etabinemc = fEtaBinInEMCalAcc;
 252   phibinemc = fPhiBinInEMCalAcc;
 253   nbinphi = fNbinPhi;
 254 }
 255
 256 //__________________________________________________________
 257 Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const
 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
 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 //__________________________________________________________
 412 void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
 413 {
 414 // Get (eta,phi) position for a given index BUT loop over all entries (takes time)
 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 //__________________________________________________________
 460 Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta)
 461 {
 462 // Get index value for a (eta,phi) position - Direct value
 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 //__________________________________________________________
 478 Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
 479 {
 480 // Get eta id
 481 // Eta discretized
 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 //__________________________________________________________
 498 Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
 499 {
 500 // Get phi id
 501 // Phi discretized
 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 //__________________________________________________________
 524 void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par)
 525 {
 526 // Allows to set parameters using only one index (if fGrid==0) !!
 527 // Not used !
 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 //__________________________________________________________
 537 void AliJetGrid::SetMatrixIndexes()
 538 {
 539 // Fill the final matrix object with the corresponding index in eta/phi
 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 //__________________________________________________________
 557 void AliJetGrid::SetIndexIJ()
 558 {
 559 // Set the grid index
 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 //__________________________________________________________
 579 void 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
 582   i = (*fIndexI)[index];
 583   j = (*fIndexJ)[index];
 584 }
 585
 586 //__________________________________________________________
 587 void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
 588 {
 589 // Returns eta, phi values for a given grid index
 590
 591   phi = fPhi->At((*fIndexI)[index]);
 592   eta = fEta->At((*fIndexJ)[index]);
 593 }
 594
 595 //__________________________________________________________
 596 Int_t AliJetGrid::GetNEntries()
 597 {
 598 // Returns the number of entries of the grid
 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 //__________________________________________________________
 612 Int_t AliJetGrid::GetNEntries2()
 613 {
 614 // Returns the number of entries of the grid
 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