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