JETAN/AliJetFillUnitArrayTracks.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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
  17 //----------------------------------------------------------------------
  18 // Fill Unit Array class
  19 // Class used by AliJetESDReader to fill a UnitArray from the information
  20 // extracted from the particle tracks
  21 // Author: magali.estienne@ires.in2p3.fr
  22 //----------------------------------------------------------------------
  23
  24
  25 // --- Standard library ---
  26 #include <Riostream.h>
  27
  28 // --- ROOT system ---
  29 #include <TSystem.h>
  30 #include <TLorentzVector.h>
  31 #include <TRefArray.h>
  32 #include <TVector3.h>
  33 #include "TTask.h"
  34 #include <TGeoManager.h>
  35 #include <TMatrixD.h>
  36 #include <TArrayD.h>
  37 #include <TMath.h>
  38 #include <TClonesArray.h>
  39
  40 // --- AliRoot header files ---
  41 #include "AliJetFinder.h"
  42 #include "AliJetReaderHeader.h"
  43 #include "AliJetReader.h"
  44 #include "AliJetESDReader.h"
  45 #include "AliJetESDReaderHeader.h"
  46 //#include "AliESD.h"
  47 #include "AliESDEvent.h"
  48 #include "AliJetDummyGeo.h"
  49 #include "AliJetUnitArray.h"
  50 #include "AliJetFillUnitArrayTracks.h"
  51 #include "AliJetHadronCorrectionv1.h"
  52 #include "AliJetGrid.h"
  53
  54 ClassImp(AliJetFillUnitArrayTracks)
  55
  56 //_____________________________________________________________________________
  57 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
  58   : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
  59     fNumUnits(0),
  60     fEtaMinCal(0),
  61     fEtaMaxCal(0),
  62     fPhiMinCal(0),
  63     fPhiMaxCal(0),
  64     fHadCorr(0),
  65     fHCorrection(0),
  66     fNTracks(0),
  67     fNTracksCut(0),
  68     fOpt(0),
  69     fDZ(0),
  70     fDebug(0),
  71     fReaderHeader(0x0),
  72     fMomentumArray(0x0),
  73     fUnitArray(0x0),
  74     fRefArray(0x0),
  75     fTPCGrid(0x0),
  76     fEMCalGrid(0x0),
  77     fGeom(0x0),
  78     fESD(0x0),
  79     fGrid0(0x0),
  80     fGrid1(0x0),
  81     fGrid2(0x0),
  82     fGrid3(0x0),
  83     fGrid4(0x0),
  84     fNphi(0),
  85     fNeta(0),
  86     fPhi2(0x0),
  87     fEta2(0x0),
  88     fPhi(0x0),
  89     fEta(0x0),
  90     fIndex(0x0),
  91     fParams(0x0),
  92     fGrid(0),
  93     fPhiMin(0),
  94     fPhiMax(0),
  95     fEtaMin(0),
  96     fEtaMax(0),
  97     fEtaBinInTPCAcc(0),
  98     fPhiBinInTPCAcc(0),
  99     fEtaBinInEMCalAcc(0),
 100     fPhiBinInEMCalAcc(0),
 101     fNbinPhi(0)
 102 {
 103   // constructor
 104 }
 105
 106 //_____________________________________________________________________________
 107 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/)
 108   : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
 109     fNumUnits(0),
 110     fEtaMinCal(0),
 111     fEtaMaxCal(0),
 112     fPhiMinCal(0),
 113     fPhiMaxCal(0),
 114     fHadCorr(0),
 115     fHCorrection(0),
 116     fNTracks(0),
 117     fNTracksCut(0),
 118     fOpt(0),
 119     fDZ(0),
 120     fDebug(0),
 121     fReaderHeader(0x0),
 122     fMomentumArray(0x0),
 123     fUnitArray(0x0),
 124     fRefArray(0x0),
 125     fTPCGrid(0x0),
 126     fEMCalGrid(0x0),
 127     fGeom(0x0),
 128     fESD(0x0),
 129     fGrid0(0x0),
 130     fGrid1(0x0),
 131     fGrid2(0x0),
 132     fGrid3(0x0),
 133     fGrid4(0x0),
 134     fNphi(0),
 135     fNeta(0),
 136     fPhi2(0x0),
 137     fEta2(0x0),
 138     fPhi(0x0),
 139     fEta(0x0),
 140     fIndex(0x0),
 141     fParams(0x0),
 142     fGrid(0),
 143     fPhiMin(0),
 144     fPhiMax(0),
 145     fEtaMin(0),
 146     fEtaMax(0),
 147     fEtaBinInTPCAcc(0),
 148     fPhiBinInTPCAcc(0),
 149     fEtaBinInEMCalAcc(0),
 150     fPhiBinInEMCalAcc(0),
 151     fNbinPhi(0)
 152 {
 153   // constructor
 154 }
 155
 156 //____________________________________________________________________________
 157 void AliJetFillUnitArrayTracks::InitParameters()
 158 {
 159   //  fHCorrection    = 0;     // For hadron correction
 160   fHadCorr        = 0;     // For hadron correction
 161   fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
 162   fDebug = fReaderHeader->GetDebug();
 163
 164   fEtaMinCal = fGeom->GetArm1EtaMin();
 165   fEtaMaxCal = fGeom->GetArm1EtaMax();
 166   fPhiMinCal = fGeom->GetArm1PhiMin();
 167   fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
 168
 169   fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
 170                         fPhiMax,fEtaMin,fEtaMax);
 171   fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
 172                         fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
 173
 174   fEta   = fTPCGrid->GetArrayEta();
 175   fPhi   = fTPCGrid->GetArrayPhi();
 176   fIndex = fTPCGrid->GetIndexObject();
 177
 178   if(fDebug>20){
 179     for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
 180     for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
 181
 182     for(Int_t i=0; i<fNphi+1; i++)
 183       for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
 184           (*fIndex)(i,j) << endl; }
 185   }
 186   if(fDebug>1) printf("\n Parameters initiated ! \n");
 187 }
 188
 189 //_____________________________________________________________________________
 190 AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
 191 {
 192   // destructor
 193 }
 194
 195 //_____________________________________________________________________________
 196 void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
 197 {
 198   //
 199   // Main method.
 200   //
 201
 202   fDebug = fReaderHeader->GetDebug();
 203
 204   // Set parameters
 205   InitParameters();
 206
 207   // get number of tracks in event (for the loop)
 208   Int_t goodTrack = 0;
 209   Int_t nt = 0;
 210   Float_t pt, eta,phi;
 211   TVector3 p3;
 212   nt = fESD->GetNumberOfTracks();
 213   if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
 214
 215  // temporary storage of signal and pt cut flag
 216   Int_t* sflag  = new Int_t[nt];
 217   Int_t* cflag  = new Int_t[nt];
 218
 219   // get cuts set by user
 220   Float_t ptMin  = fReaderHeader->GetPtCut();
 221   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 222   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 223   fOpt = fReaderHeader->GetDetector();
 224   fDZ  = fReaderHeader->GetDZ();
 225
 226   Int_t nTracksEmcal      = 0;
 227   Int_t nTracksEmcalDZ    = 0;
 228   Int_t nTracksTpc        = 0;
 229   Int_t nTracksTpcOnly    = 0;
 230   Int_t nTracksEmcalCut   = 0;
 231   Int_t nTracksEmcalDZCut = 0;
 232   Int_t nTracksTpcCut     = 0;
 233   Int_t nTracksTpcOnlyCut = 0;
 234
 235   fGrid = fTPCGrid->GetGridType();
 236
 237
 238   //loop over tracks
 239   for (Int_t it = 0; it < nt; it++) {
 240     AliESDtrack *track = fESD->GetTrack(it);
 241     UInt_t status = track->GetStatus();
 242
 243     Double_t mom[3];
 244     track->GetPxPyPz(mom);
 245     p3.SetXYZ(mom[0],mom[1],mom[2]);
 246     pt = p3.Pt();
 247     Float_t mass = 0.;
 248     mass = track->GetMass();
 249
 250     if (((status & AliESDtrack::kITSrefit) == 0) ||
 251         ((status & AliESDtrack::kTPCrefit) == 0)) continue;    // quality check
 252     if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
 253         && TMath::Abs(track->GetLabel()) > 10000)  continue;   // quality check
 254     if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
 255         && TMath::Abs(track->GetLabel()) < 10000)  continue;   // quality check
 256     eta = p3.Eta();
 257     phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
 258
 259     if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
 260
 261     // sflag -> not yet implemented !!!!
 262
 263     if(fGrid==0)
 264       {
 265         // Only TPC filled from its grid in its total acceptance
 266
 267         Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
 268         Bool_t ok = kFALSE;
 269
 270         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
 271         uArray->SetUnitTrackID(it);
 272
 273         Float_t unitEnergy = 0.;
 274         unitEnergy = uArray->GetUnitEnergy();
 275         if(unitEnergy==0.){
 276           nTracksTpcOnly++;
 277           ok = kTRUE;
 278         }
 279         // Fill energy in TPC acceptance
 280         uArray->SetUnitEnergy(unitEnergy + pt);
 281         uArray->SetUnitPxPyPz(mom);
 282         uArray->SetUnitMass(mass);
 283
 284         // Pt cut flag
 285         if(uArray->GetUnitEnergy()<ptMin){
 286           uArray->SetUnitCutFlag(kPtSmaller);
 287         }
 288         else {
 289           uArray->SetUnitCutFlag(kPtHigher);
 290           if(ok) nTracksTpcOnlyCut++;
 291         }
 292
 293         // Detector flag
 294         if(unitEnergy>0) {
 295           uArray->SetUnitDetectorFlag(kAll);
 296         }
 297         if(uArray->GetUnitEnergy()>0){
 298           fRefArray->Add(uArray);
 299         }
 300
 301         sflag[goodTrack]=0;
 302         if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 303         cflag[goodTrack]=0;
 304         if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
 305         goodTrack++;
 306
 307       }
 308
 309     if(fGrid==1)
 310       {
 311         Int_t nElements = fTPCGrid->GetNEntries();
 312
 313         // Fill track information in EMCAL acceptance
 314         if((eta >= fEtaMin && eta <= fEtaMax) &&
 315            (phi >= fPhiMin && phi <= fPhiMax))// &&
 316           {
 317
 318             // Include dead-zones
 319             if(fDZ)
 320               {
 321                 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 322                 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 323                 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 324                 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 325                 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 326                 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 327                 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 328                 Int_t n0 = fGrid0->GetNEntries();
 329                 Int_t n1 = fGrid1->GetNEntries();
 330                 Int_t n2 = fGrid2->GetNEntries();
 331                 Int_t n3 = fGrid3->GetNEntries();
 332
 333                 if(phi >= phimin0 && phi <= phimax0){
 334                   Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
 335                   AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
 336                   uArray0->SetUnitTrackID(it);
 337                   Float_t uEnergy0 = uArray0->GetUnitEnergy();
 338                   Bool_t ok0 = kFALSE;
 339                   if(uEnergy0==0.){
 340                     nTracksEmcalDZ++;
 341                     ok0 = kTRUE;
 342                   }
 343                   uArray0->SetUnitEnergy(uEnergy0+pt);
 344                   if(uArray0->GetUnitEnergy()<ptMin)
 345                     uArray0->SetUnitCutFlag(kPtSmaller);
 346                   else {
 347                     uArray0->SetUnitCutFlag(kPtHigher);
 348                     if(ok0) nTracksEmcalDZCut++;
 349                   }
 350                   if(uArray0->GetUnitEnergy()>0)
 351                     fRefArray->Add(uArray0);
 352                 }
 353                 if(phi >= phimin1 && phi <= phimax1){
 354                   Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
 355                   AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
 356                   uArray1->SetUnitTrackID(it);
 357                   Float_t uEnergy1 = uArray1->GetUnitEnergy();
 358                   Bool_t ok1 = kFALSE;
 359                   if(uEnergy1==0.){
 360                     nTracksEmcalDZ++;
 361                     ok1 = kTRUE;
 362                   }
 363                   uArray1->SetUnitEnergy(uEnergy1+pt);
 364                   if(uArray1->GetUnitEnergy()<ptMin)
 365                     uArray1->SetUnitCutFlag(kPtSmaller);
 366                   else {
 367                     uArray1->SetUnitCutFlag(kPtHigher);
 368                     if(ok1) nTracksEmcalDZCut++;
 369                   }
 370                   if(uArray1->GetUnitEnergy()>0)
 371                     fRefArray->Add(uArray1);
 372                 }
 373                 if(phi >= phimin2 && phi <= phimax2){
 374                   Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
 375                   AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
 376                   uArray2->SetUnitTrackID(it);
 377                   Float_t uEnergy2 = uArray2->GetUnitEnergy();
 378                   Bool_t ok2 = kFALSE;
 379                   if(uEnergy2==0.){
 380                     nTracksEmcalDZ++;
 381                     ok2 = kTRUE;
 382                   }
 383                   uArray2->SetUnitEnergy(uEnergy2+pt);
 384                   if(uArray2->GetUnitEnergy()<ptMin)
 385                     uArray2->SetUnitCutFlag(kPtSmaller);
 386                   else {
 387                     uArray2->SetUnitCutFlag(kPtHigher);
 388                     if(ok2) nTracksEmcalDZCut++;
 389                   }
 390                   if(uArray2->GetUnitEnergy()>0)
 391                     fRefArray->Add(uArray2);
 392                 }
 393                 if(phi >= phimin3 && phi <= phimax3){
 394                   Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
 395                   AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
 396                   uArray3->SetUnitTrackID(it);
 397                   Float_t uEnergy3 = uArray3->GetUnitEnergy();
 398                   Bool_t ok3 = kFALSE;
 399                   if(uEnergy3==0.){
 400                     nTracksEmcalDZ++;
 401                     ok3 = kTRUE;
 402                   }
 403                   uArray3->SetUnitEnergy(uEnergy3+pt);
 404                   if(uArray3->GetUnitEnergy()<ptMin)
 405                     uArray3->SetUnitCutFlag(kPtSmaller);
 406                   else {
 407                     uArray3->SetUnitCutFlag(kPtHigher);
 408                     if(ok3) nTracksEmcalDZCut++;
 409                   }
 410                   if(uArray3->GetUnitEnergy()>0)
 411                     fRefArray->Add(uArray3);
 412                 }
 413                 if(phi >= phimin4 && phi <= phimax4){
 414                   Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
 415                   AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
 416                   uArray4->SetUnitTrackID(it);
 417                   Float_t uEnergy4 = uArray4->GetUnitEnergy();
 418                   Bool_t ok4 = kFALSE;
 419                   if(uEnergy4==0.){
 420                     nTracksEmcalDZ++;
 421                     ok4 = kTRUE;
 422                   }
 423                   uArray4->SetUnitEnergy(uEnergy4+pt);
 424                   if(uArray4->GetUnitEnergy()<ptMin)
 425                     uArray4->SetUnitCutFlag(kPtSmaller);
 426                   else {
 427                     uArray4->SetUnitCutFlag(kPtHigher);
 428                     if(ok4) nTracksEmcalDZCut++;
 429                   }
 430                   if(uArray4->GetUnitEnergy()>0)
 431                     fRefArray->Add(uArray4);
 432                 }
 433               } // end fDZ
 434
 435             Int_t towerID = 0;
 436             fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
 437
 438             if(towerID==-1) continue;
 439
 440             AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
 441             uArray->SetUnitTrackID(it);
 442             Float_t unitEnergy = uArray->GetUnitEnergy();
 443             Bool_t ok = kFALSE;
 444             if(unitEnergy==0.){
 445               nTracksEmcal++;
 446               ok=kTRUE;
 447             }
 448
 449             // Do Hadron Correction
 450             // Parametrization to be added
 451             if (fHCorrection != 0)
 452               {
 453                 //            Float_t   hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
 454                 Float_t   hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta);
 455                 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
 456
 457               } //end Hadron Correction loop
 458
 459             uArray->SetUnitEnergy(unitEnergy + pt);
 460
 461             // Put a pt cut flag
 462             if(uArray->GetUnitEnergy()<ptMin)
 463               uArray->SetUnitCutFlag(kPtSmaller);
 464             else {
 465               uArray->SetUnitCutFlag(kPtHigher);
 466               if(ok) nTracksEmcalCut++;
 467             }
 468             // Detector flag
 469             if(unitEnergy > 0)
 470               uArray->SetUnitDetectorFlag(kAll);
 471
 472             if(uArray->GetUnitEnergy()>0)
 473               fRefArray->Add(uArray);
 474
 475             sflag[goodTrack]=0;
 476             if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 477             cflag[goodTrack]=0;
 478             if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
 479             goodTrack++;
 480
 481           } // end loop on EMCal acceptance cut + tracks quality
 482         else{
 483           // Outside EMCal acceptance
 484
 485           //    Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
 486           Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
 487
 488           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
 489           uArray->SetUnitTrackID(it);
 490
 491           Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
 492           Bool_t ok2 = kFALSE;
 493           if(unitEnergy2==0.){
 494             nTracksTpc++;
 495             ok2=kTRUE;
 496           }
 497           // Fill energy outside emcal acceptance
 498           uArray->SetUnitEnergy(unitEnergy2 + pt);
 499
 500           // Pt cut flag
 501           if(uArray->GetUnitEnergy()<ptMin){
 502             uArray->SetUnitCutFlag(kPtSmaller);
 503           }
 504           else {
 505             uArray->SetUnitCutFlag(kPtHigher);
 506             if(ok2) nTracksTpcCut++;
 507           }
 508           // Detector flag
 509           if(unitEnergy2 > 0)
 510             uArray->SetUnitDetectorFlag(kTpc);
 511           if(uArray->GetUnitEnergy()>0)
 512             fRefArray->Add(uArray);
 513
 514           sflag[goodTrack]=0;
 515           if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 516           cflag[goodTrack]=0;
 517           if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
 518           goodTrack++;
 519
 520         }
 521       } // end fGrid==1
 522   } // end loop on entries (tpc tracks)
 523
 524 //   // set the signal flags
 525 //   fSignalFlag.Set(goodTrack,sflag);
 526 //   fCutFlag.Set(goodTrack,cflag);
 527
 528 //   delete sflag;
 529 //   delete cflag;
 530
 531   //    } // end loop on entries (tpc tracks)
 532
 533   if(fGrid==0) {
 534     fNTracks = nTracksTpcOnly;
 535     fNTracksCut = nTracksTpcOnlyCut;
 536     if(fDebug>10){
 537       cout << "fNTracks : " << fNTracks << endl;
 538       cout << "fNTracksCut : " << fNTracksCut << endl;
 539     }
 540   }
 541   if(fGrid==1) {
 542     fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
 543     fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
 544     if(fDebug>10){
 545       cout << "fNTracks : " << fNTracks << endl;
 546       cout << "fNTracksCut : " << fNTracksCut << endl;
 547     }
 548   }
 549
 550 }
 551
 552 //__________________________________________________________
 553 void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
 554 {
 555 //
 556 // Obtain the index from eta and phi
 557 //
 558   for(Int_t j=0; j<fNphi+1; j++) {
 559     for(Int_t i=0; i<fNeta+1; i++) {
 560
 561       // TPC grid only
 562       //-------------------------------------
 563       if(fGrid==0) {
 564         if(j*(fNeta+1)+i == index) {
 565           eta = fEta2->At(i);
 566           phi = fPhi2->At(j);
 567         }
 568       }
 569
 570       // TPC-EMCAL grid
 571       //-------------------------------------
 572       Int_t ii = 0;
 573       if(i==0) ii = 0;
 574       if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
 575       if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
 576
 577       if(fGrid==1) {
 578         if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
 579           eta = fEta2->At(i);
 580           phi = fPhi2->At(j);
 581         }
 582
 583         if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
 584            ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
 585           if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
 586           else eta = fEta2->At(i);
 587           phi = fPhi2->At(j);
 588         }
 589
 590         if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
 591           eta = fEta2->At(i);
 592           phi = fPhi2->At(j);
 593         }
 594       }
 595     }
 596   }
 597 }
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609