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 { // To set initial parameters
 168
 169   fPhiMin = phiMinCal; // rad
 170   fPhiMax = phiMaxCal; // rad
 171   fEtaMin = etaMinCal;
 172   fEtaMax = etaMaxCal;
 173   fNbinPhi = static_cast<int>(fPhiMin/(2*TMath::Pi()/fNphi));
 174
 175   // Define some binning numbers
 176   if(fGrid==0){
 177     for(Int_t i=0; i<fNphi+1; i++) fPhiBinInTPCAcc++;
 178     fPhiBinInEMCalAcc = 0;
 179
 180     for(Int_t i=0; i<fNeta+1; i++) fEtaBinInTPCAcc++;
 181     fEtaBinInEMCalAcc = 0;
 182   }
 183
 184   if(fGrid==1){
 185     for(Int_t i=0; i<fNphi+1; i++) {
 186       fPhiBinInTPCAcc++;
 187       if(fPhi->At(i) >= fPhiMin &&
 188          fPhi->At(i) <= fPhiMax)
 189         fPhiBinInEMCalAcc++;
 190     }
 191     for(Int_t i=0; i<fNeta+1; i++) {
 192       fEtaBinInTPCAcc++;
 193       if((fEta->At(i) >= fEtaMin) &&
 194          (fEta->At(i) <= fEtaMax))
 195         fEtaBinInEMCalAcc++;
 196     }
 197   }
 198
 199 }
 200
 201 //__________________________________________________________
 202 TArrayD* AliJetGrid::GetArrayEta()
 203 { // Returns an array with the eta points
 204
 205   return fEta;
 206
 207 }
 208
 209 //__________________________________________________________
 210 TArrayD* AliJetGrid::GetArrayPhi()
 211 { // Returns an array with the phi points
 212
 213   return fPhi;
 214
 215 }
 216
 217 //__________________________________________________________
 218 TMatrixD* AliJetGrid::GetIndexObject()
 219 { // Returns a pointer to the matrix
 220
 221   return fIndex;
 222
 223 }
 224
 225 //__________________________________________________________
 226 void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi,
 227                                 Float_t &mineta, Float_t &maxeta)
 228 { // Returns EMCAL acceptance initially setted
 229
 230   nphi = fNphi;
 231   neta = fNeta;
 232   minphi = fPhiMin;
 233   maxphi = fPhiMax;
 234   mineta = fEtaMin;
 235   maxeta = fEtaMax;
 236 }
 237
 238 //__________________________________________________________
 239 void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc,
 240                                 Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi)
 241 { // Returns number of bins in TPC and EMCAL geometry
 242
 243   etabintpc = fEtaBinInTPCAcc;
 244   phibintpc = fPhiBinInTPCAcc;
 245   etabinemc = fEtaBinInEMCalAcc;
 246   phibinemc = fPhiBinInEMCalAcc;
 247   nbinphi = fNbinPhi;
 248 }
 249
 250 //__________________________________________________________
 251 Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const
 252 { // Tells the index value of a corresponding (eta,phi) real position
 253   // Loop over all entries -> takes time.
 254   // Used one time at the begining to fill the grids
 255
 256   /*   this is how bins are numbered
 257
 258        in all TPC        or      in TPC - EMCAL
 259                                 ... ... ... ..
 260                                 ---------------
 261           ...  ... .            10 |   |   | 11
 262           ---+---+---           ---------------
 263     ^      6 | 7 | 8     or      8 |   |   | 9
 264     |     ---+---+---           ---------------
 265            3 | 4 | 5             4 | 5 | 6 | 7
 266    phi    ---+---+---           ---------------
 267            0 | 1 | 2             0 | 1 | 2 | 3
 268
 269
 270              eta ->
 271   */
 272
 273   Int_t etaBin=0,phiBin=0,absID=0;
 274   Int_t etaBin2=0,etaBin3=0;
 275
 276   // Fill all the grid in eta/phi (all TPC acceptance)
 277   //-----------------------------------------------------
 278   if(fGrid==0){
 279     if(eta <= fEta->At(0)) {
 280       etaBin = 0;
 281     } else if(eta >= fEta->At(fNeta)) {
 282       etaBin = fNeta;
 283     } else {
 284       for(Int_t i=0; i<fNeta+1; i++) {
 285         if(eta < fEta->At(i)) {
 286           etaBin = i-1;
 287           break;
 288         }
 289       }
 290     }
 291     if(phi <= fPhi->At(0)) {
 292       phiBin = 0;
 293     } else if(phi >= fPhi->At(fNphi)) {
 294       phiBin = fNphi;
 295     } else {
 296       for(Int_t i=0; i<fNphi+1; i++) {
 297         if(phi < fPhi->At(i)) {
 298           phiBin = i-1;
 299           break;
 300         }
 301       }
 302     }
 303
 304     // Calculate absolute id
 305     absID = phiBin*(fNeta+1) + etaBin;
 306
 307   }
 308
 309   // Fill the grid but do not count id in EMCal acceptance
 310   //------------------------------------------------------
 311   if((eta >= fEtaMin && eta <= fEtaMax) &&
 312      (phi >= fPhiMin && phi <= fPhiMax)){
 313     etaBin = etaBin2 = etaBin3 = -1;
 314     phiBin = -1;
 315   }
 316
 317   if(fGrid == 1){
 318     if(phi<fPhiMin) {
 319       if(eta <= fEta->At(0)) {
 320         etaBin = 0;
 321       } else if(eta >= fEta->At(fNeta)) {
 322         etaBin = fNeta;
 323       } else {
 324         for(Int_t i=0; i<fNeta+1; i++) {
 325           if(eta < fEta->At(i)) {
 326             etaBin = i-1;
 327             break;
 328           }
 329         }
 330       }
 331     }
 332     else {
 333       if((phi>=fPhiMin) && (phi<=fPhiMax)) {
 334         if(eta <= fEta->At(0)) {
 335           etaBin2 = 0;
 336         } else if(eta >= fEta->At(fNeta)) {
 337           etaBin2 = fNeta-fEtaBinInEMCalAcc;
 338         } else {
 339           for(Int_t i=0; i<fNeta+1; i++) {
 340             if(eta < fEta->At(i) && ((fEta->At(0)<eta && eta<fEtaMin) ||
 341                                      (fEtaMax<eta && eta<fEta->At(fNeta)))) {
 342               //                 cout << "i : " << i << endl;
 343               if(eta<fEtaMin)
 344                 etaBin2 = i-1;
 345               else etaBin2 = i-1-fEtaBinInEMCalAcc;
 346               break;
 347             }
 348           }
 349         }
 350       }
 351       else {
 352         if(eta <= fEta->At(0)) {
 353           etaBin3 = 0;
 354         } else if(eta >= fEta->At(fNeta)) {
 355           etaBin3 = fNeta;
 356         } else {
 357           for(Int_t i=0; i<fNeta+1; i++) {
 358             if(eta < fEta->At(i)) {
 359               etaBin3 = i-1;
 360               break;
 361             }
 362           }
 363         }
 364       }
 365     }
 366
 367     if(phi <= fPhi->At(0)) {
 368       phiBin = 0;
 369     } else if(phi >= fPhi->At(fNphi)) {
 370       phiBin = fNphi;
 371     } else {
 372       for(Int_t i=0; i<fNphi+1; i++) {
 373         if(phi < fPhi->At(i)) {
 374           phiBin = i-1;
 375           break;
 376         }
 377       }
 378     }
 379
 380     // Calculate absID
 381     //-------------------------------------------------------
 382     if(phi<fPhiMin)
 383       absID = phiBin*(fNeta+1) + etaBin;
 384     if(phi>=fPhiMin && phi<=fPhiMax) {
 385       if(eta >= fEtaMin && eta <= fEtaMax) absID = -1;
 386       else{
 387         absID = (fNbinPhi+1)*(fNeta+1)
 388           + (phiBin-fNbinPhi-1)*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
 389           + etaBin2;
 390       }
 391     }
 392     if(phi>fPhiMax)
 393       absID = (fNbinPhi+1)*(fNeta+1)
 394         + fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
 395         + (phiBin-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)
 396         + etaBin3;
 397
 398   } // END OPTION==1
 399
 400   return absID;
 401
 402 }
 403
 404 //__________________________________________________________
 405 void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
 406 { // Get (eta,phi) position for a given index BUT loop over all entries (takes time)
 407
 408   for(Int_t j=0; j<fNphi+1; j++) {
 409     for(Int_t i=0; i<fNeta+1; i++) {
 410
 411       // TPC grid only
 412       //-------------------------------------
 413       if(fGrid==0) {
 414         if(j*(fNeta+1)+i == index) {
 415           eta = fEta->At(i);
 416           phi = fPhi->At(j);
 417         }
 418       }
 419
 420       // TPC-EMCAL grid
 421       //-------------------------------------
 422       Int_t ii = 0;
 423       if(i==0) ii = 0;
 424       if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
 425       if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
 426
 427       if(fGrid==1) {
 428         if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
 429           eta = fEta->At(i);
 430           phi = fPhi->At(j);
 431         }
 432
 433         if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
 434            ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
 435           if(ii==0) {Int_t ind = 0; eta = fEta->At(ind);}
 436           else eta = fEta->At(i);
 437           phi = fPhi->At(j);
 438         }
 439
 440         if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) &&
 441            ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
 442             +(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
 443           eta = fEta->At(i);
 444           phi = fPhi->At(j);
 445         }
 446       }
 447     }
 448   }
 449 }
 450
 451 //__________________________________________________________
 452 Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta)
 453 { // Get index value for a (eta,phi) position - Direct value
 454
 455   Int_t ieta = GetIndexJFromEta(eta);
 456   Int_t iphi = GetIndexIFromPhi(phi);
 457
 458   Int_t index = GetMatrixIndex(iphi,ieta);
 459
 460   if(fDebug>10){
 461     cout << "(phi,eta) : " << phi << ", " << eta << endl;
 462     cout << "index : " << index << endl;
 463   }
 464   return index;
 465
 466 }
 467
 468 //__________________________________________________________
 469 Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
 470 { // Get eta id
 471   // Eta discretized
 472
 473   Int_t idEta =0;
 474   Double_t temp = (eta+fMaxEta)/(fMaxEta-fMinEta)*fNeta;
 475   if(fDebug>20)
 476     {
 477       cout << "eta : " << eta << endl;
 478       cout << "fMaxEta : " << fMaxEta << endl;
 479       cout << "fMinEta : " << fMinEta << endl;
 480       cout << "fNeta : " << fNeta << endl;
 481       cout << "temp eta before cast : " << temp << endl;
 482     }
 483   idEta = static_cast<Int_t>(temp+0.5);
 484   if(fDebug>20) cout << "temp eta after cast : " << idEta << endl;
 485   return idEta;
 486 }
 487 //__________________________________________________________
 488 Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
 489 { // Get phi id
 490   // Phi discretized
 491
 492   Int_t idPhi = 0;
 493   Double_t temp = 0.;
 494   if(fMinPhi==0) temp = phi/(fMaxPhi-fMinPhi)*fNphi;
 495   else temp = (phi-fMinPhi)/(fMaxPhi-fMinPhi)*fNphi;
 496
 497   if(fDebug>20)
 498     {
 499       cout << "phi : " << phi << endl;
 500       cout << "fMaxPhi : " << fMaxPhi << endl;
 501       cout << "fMinPhi : " << fMinPhi << endl;
 502       cout << "fNphi : " << fNphi << endl;
 503       cout << "temp phi before cast : " << temp << endl;
 504     }
 505   idPhi = static_cast<Int_t>(temp+0.5);
 506   if(fDebug>20) cout << "temp phi after cast : " << idPhi << endl;
 507   return idPhi;
 508
 509
 510 }
 511
 512 //__________________________________________________________
 513 void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par)
 514 { // Allows to set parameters using only one index (if fGrid==0) !!
 515   // Not used !
 516
 517   Int_t iphi = (Int_t)i/fNeta;
 518   Int_t ieta = i-iphi*fNeta;
 519   SetMatrixIndex(iphi,ieta,par);
 520
 521   return;
 522 }
 523
 524 //__________________________________________________________
 525 void AliJetGrid::SetMatrixIndexes()
 526 { // Fill the final matrix object with the corresponding index in eta/phi
 527
 528   for(Int_t i=0; i<fNphi+1; i++){
 529     for(Int_t j=0; j<fNeta+1; j++){
 530       (*fIndex)(i,j) = GetIndexFromEtaPhi(fPhi->At(i),fEta->At(j))+1;
 531       if(fDebug>2){
 532         cout << "(*fIndex)(" << i << "," << j << ") : " << (*fIndex)(i,j) <<
 533           ", phi : " << fPhi->At(i) << ", eta : " << fEta->At(j) << endl;
 534       }
 535     }
 536   }
 537   printf("##########################################################\n");
 538   printf("TMatrix object filled !\n");
 539   printf("Size of the object : phi x eta = (fNphi+1) x (fNeta+1) = %d\n",(fNphi+1)*(fNeta+1));
 540   printf("##########################################################\n");
 541 }
 542
 543 //__________________________________________________________
 544 void AliJetGrid::SetIndexIJ()
 545 { //
 546
 547   for(Int_t i=0; i<fNphi+1; i++){
 548     for(Int_t j=0; j<fNeta+1; j++){
 549       Int_t id = static_cast<Int_t>((*fIndex)(i,j));
 550
 551       if(id!=-1)
 552         {
 553           (*fIndexI)[id] = i;
 554           (*fIndexJ)[id] = j;
 555         }
 556     }
 557   }
 558
 559   printf("##########################################################\n");
 560   printf("     In SetIndexIJ - Grid indexes setted !\n");
 561   printf("##########################################################\n");
 562 }
 563
 564 //__________________________________________________________
 565 void AliJetGrid::GetIJFromIndex(Int_t index, Int_t i, Int_t j)
 566 { // returns i position id of eta and j position id of phi for a given grid index
 567   i = (*fIndexI)[index];
 568   j = (*fIndexJ)[index];
 569 }
 570
 571 //__________________________________________________________
 572 void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
 573 { // returns eta, phi values for a given grid index
 574
 575   phi = fPhi->At((*fIndexI)[index]);
 576   eta = fEta->At((*fIndexJ)[index]);
 577 }
 578
 579 //__________________________________________________________
 580 Int_t AliJetGrid::GetNEntries()
 581 { // Returns the number of entries of the grid
 582
 583   if(fDebug>20){
 584     cout << "fMaxPhi : " << fMaxPhi << endl;
 585     cout << "fMaxEta : " << fMaxEta << endl;
 586   }
 587
 588   Int_t indexNum = GetIndex(fMaxPhi,fMaxEta);
 589   if(fDebug>20) cout << "indexNum : " << indexNum << endl;
 590   return indexNum;
 591
 592 }
 593
 594 //__________________________________________________________
 595 Int_t AliJetGrid::GetNEntries2()
 596 { // Returns the number of entries of the grid
 597
 598   Int_t indexNum = GetIndex(fMaxPhi-1.,fMaxEta-0.5);
 599     if(fDebug>20) cout << "indexNum : " << indexNum << endl;
 600   return indexNum;
 601
 602 }
 603
 604
 605
 606
 607
 608