JETAN/AliJetFillUnitArrayTracks.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //======================================================================
  18 // *** November 2006
  19 // Author: magali.estienne@ires.in2p3.fr
  20 // 1) Define 2 grids and take (eta,phi) from the grid or use the grid for the TPC and
  21 // EtaPhiFromIndex and TowerIndexFromEtaPhi for the particles in EMCAL acceptance
  22 //     2 options are possible : fGrid==0, work only with full TPC acceptance (for the moment)
  23 //                              fGrid==1, work with a part of the TPC acceptance
  24 // 2) Try to implement 2 full grids for TPC and EMCal separately and to merge them
  25 // 3) Need to include Dead-zone -> Wait for exact positions in the new detector geometry
  26 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
  27 //======================================================================
  28 // ***September 2006
  29 // TTask : Fill Unit Array for the Tracks information
  30 // Called by ESD reader for jet analysis
  31 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
  32 //======================================================================
  33 // *** July 2006
  34 // 1) When the tracks are in the EMCal acceptance, the functions EtaPhiFromIndex
  35 // and TowerIndexFromEtaPhi in the AliEMCALGeometry class are used to extract the
  36 // index or the eta, phi position of a grid.
  37 // 2) Define a grid for TPC
  38 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
  39 //======================================================================
  40 // ***July 2006
  41 // Fill Unit Array class
  42 // Class used by AliJetESDReader to fill a UnitArray from the information extracted
  43 // from the particle tracks
  44 // Author: magali.estienne@ires.in2p3.fr
  45 //======================================================================
  46
  47
  48 // --- Standard library ---
  49 #include <Riostream.h>
  50
  51 // --- ROOT system ---
  52 #include <TSystem.h>
  53 #include <TLorentzVector.h>
  54 #include <TVector3.h>
  55 //#include "Math/Vector3D.h"
  56 //#include "Math/Vector3Dfwd.h"
  57 #include "TTask.h"
  58 #include <TGeoManager.h>
  59 #include <TMatrixD.h>
  60 #include <TArrayD.h>
  61
  62 // --- AliRoot header files ---
  63 #include "AliJetFinder.h"
  64 #include "AliJetReaderHeader.h"
  65 #include "AliJetReader.h"
  66 #include "AliJetESDReader.h"
  67 #include "AliJetESDReaderHeader.h"
  68 #include "AliESD.h"
  69 //#include "AliEMCALGeometry.h"
  70 #include "AliJetDummyGeo.h"
  71 #include "AliJetUnitArray.h"
  72 #include "AliJetFillUnitArrayTracks.h"
  73 #include "AliJetHadronCorrectionv1.h"
  74 #include "AliJetGrid.h"
  75
  76 ClassImp(AliJetFillUnitArrayTracks)
  77
  78 //_____________________________________________________________________________
  79 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
  80     : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
  81       fNumUnits(0),
  82       fEtaMinCal(0.),
  83       fEtaMaxCal(0.),
  84       fPhiMinCal(0.),
  85       fPhiMaxCal(0.),
  86       fHadCorr(0x0),
  87       fHCorrection(0),
  88       fNIn(0),
  89       fOpt(0),
  90       fDebug(0),
  91       fReaderHeader(0x0),
  92       fMomentumArray(0x0),
  93       fUnitArray(0x0),
  94       fTPCGrid(0x0),
  95       fEMCalGrid(0x0),
  96       fGeom(0x0),
  97       fNphi(0),
  98       fNeta(0),
  99       fPhi2(0x0),
 100       fEta2(0x0),
 101       fPhi(0x0),
 102       fEta(0x0),
 103       fIndex(0x0),
 104       fParams(0x0),
 105       fGrid(0),
 106       fPhiMin(0.),
 107       fPhiMax(0.),
 108       fEtaMin(0.),
 109       fEtaMax(0.),
 110       fEtaBinInTPCAcc(0),
 111       fPhiBinInTPCAcc(0),
 112       fEtaBinInEMCalAcc(0),
 113       fPhiBinInEMCalAcc(0),
 114       fNbinPhi(0)
 115 {
 116   // constructor
 117 }
 118
 119 //____________________________________________________________________________
 120 void AliJetFillUnitArrayTracks::SetEMCALGeometry()
 121 {
 122     // Set EMCAL geometry information
 123     fGeom = AliJetDummyGeo::GetInstance();
 124     if (fGeom == 0)
 125         fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
 126     if(fDebug>1) printf("\n EMCAL Geometry setted ! \n");
 127 }
 128
 129 //____________________________________________________________________________
 130 void AliJetFillUnitArrayTracks::InitParameters()
 131 {
 132   fHCorrection    = 0;     // For hadron correction
 133   fHadCorr        = 0;     // For hadron correction
 134   fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
 135   fDebug = fReaderHeader->GetDebug();
 136
 137   fEtaMinCal = fGeom->GetArm1EtaMin();
 138   fEtaMaxCal = fGeom->GetArm1EtaMax();
 139   fPhiMinCal = fGeom->GetArm1PhiMin();
 140   fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
 141
 142   if(fDebug>30){
 143     cout << "fEtaMinCal : " << fEtaMinCal << endl;
 144     cout << "fEtaMaxCal : " << fEtaMaxCal << endl;
 145     cout << "fPhiMinCal : " << fPhiMinCal << endl;
 146     cout << "fPhiMaxCal : " << fPhiMaxCal << endl;
 147   }
 148
 149   fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
 150                         fPhiMax,fEtaMin,fEtaMax);
 151   fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
 152                         fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
 153
 154   fEta   = fTPCGrid->GetArrayEta();
 155   fPhi   = fTPCGrid->GetArrayPhi();
 156   fIndex = fTPCGrid->GetIndexObject();
 157
 158   if(fDebug>20){
 159     for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
 160     for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
 161
 162     for(Int_t i=0; i<fNphi+1; i++)
 163       for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
 164           (*fIndex)(i,j) << endl; }
 165   }
 166   if(fDebug>1) printf("\n Parameters initiated ! \n");
 167 }
 168
 169 //_____________________________________________________________________________
 170 AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
 171 {
 172   // destructor
 173 }
 174
 175 //_____________________________________________________________________________
 176 void AliJetFillUnitArrayTracks::Exec(Option_t* option)
 177 {
 178   //
 179   // Main method.
 180   // Explain
 181
 182   fDebug = fReaderHeader->GetDebug();
 183   if(fDebug>1) printf("In AliJetFillUnitArrayTracks !");
 184   if(fDebug>3) printf("\nfGeom->GetEntries() = %d\n", fGeom->GetNCells());
 185   // Set EMCal Geometry
 186   SetEMCALGeometry();
 187   // Set parameters
 188   InitParameters();
 189
 190   TClonesArray *lvArray = fMomentumArray;     // Correct checked !
 191   Int_t nInT =  lvArray->GetEntries();        // Correct checked !
 192   Float_t ptMin = fReaderHeader->GetPtCut();  // Correct checked !
 193
 194   // sflag -> not yet implemented !!!!
 195
 196   if(fDebug>3) cout << "nInT : " << nInT << endl;
 197
 198   if (nInT == 0) return;
 199
 200   // local arrays for input
 201   Float_t* enT  = new Float_t[nInT];
 202   Float_t* ptT  = new Float_t[nInT];
 203   Float_t* etaT = new Float_t[nInT];
 204   Float_t* phiT = new Float_t[nInT];
 205   Float_t* thetaT = new Float_t[nInT];
 206
 207   Int_t trackInEmcalAcc = 0;
 208   Int_t trackInTpcAcc = 0;
 209   Int_t trackInTpcAccOnly = 0;
 210
 211   Int_t nElements = fTPCGrid->GetNEntries();
 212   Int_t nElements2 = fEMCalGrid->GetNEntries();
 213   fGrid = fTPCGrid->GetGridType();
 214
 215   if(fDebug>3){
 216     cout << "nElements : " << nElements << endl;
 217     cout << "nElements2 : " << nElements2 << endl;
 218     cout << "fNumUnits : " << fNumUnits << endl;
 219     cout << "sum : " << fNumUnits+nElements << endl;
 220   }
 221
 222   // Set energy exactly to zero
 223   if(fGrid==0)
 224     for(Int_t k=0; k<nElements; k++)
 225       fUnitArray[k].SetUnitEnergy(0.);
 226
 227   if(fGrid==1)
 228     for(Int_t k=0; k<fNumUnits+nElements; k++)
 229       fUnitArray[k].SetUnitEnergy(0.);
 230
 231   // load input vectors
 232   for (Int_t i = 0; i < nInT; i++)
 233     {
 234       TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
 235       enT[i]  = lv->Energy();
 236       ptT[i]  = lv->Pt();
 237       etaT[i] = lv->Eta();
 238       phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
 239       thetaT[i] = 2.0*TMath::ATan(TMath::Exp(-etaT[i]));
 240
 241       if(fDebug>20){
 242         cout << "enT[" << i << "] : " <<  enT[i] << endl;
 243         cout << "ptT[" << i << "] : " <<  ptT[i] << endl;
 244         cout << "etaT[" << i << "] : " <<  etaT[i] << endl;
 245         cout << "phiT[" << i << "] : " <<  phiT[i] << endl;
 246         cout << "thetaT[" << i << "] : " <<  thetaT[i] << endl;
 247         cout << "fEtaMinCal : " << fEtaMinCal << ", fEtaMaxCal : " << fEtaMaxCal << endl;
 248         cout << "fPhiMinCal : " << fPhiMinCal << ", fPhiMaxCal : " << fPhiMaxCal << endl;
 249         cout << "fEtaMin : " << fEtaMin << ", fEtaMax : " << fEtaMax << endl;
 250         cout << "fPhiMin : " << fPhiMin << ", fPhiMax : " << fPhiMax << endl;
 251       }
 252
 253       if(fGrid==0)
 254         {
 255           // For the moment, only TPC filled from its grid in its total acceptance
 256           if(fDebug>2)
 257             cout << "In total TPC acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
 258
 259             trackInTpcAccOnly += 1;
 260
 261             Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
 262
 263             Float_t unitEnergy = 0.;
 264             unitEnergy = fUnitArray[idTPC].GetUnitEnergy();
 265
 266             if(unitEnergy > 0. && fDebug >10){
 267               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 268               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 269               cout << "TPC id : " << idTPC << endl;
 270               cout << "etaT[" << i << "] : " << etaT[i] << endl;
 271               cout << "phiT[" << i << "] : " << phiT[i] << endl;
 272               cout << "unitEnergy in TPC acceptance : " << unitEnergy << endl;
 273               cout << "fUnitArray[idTPC].GetUnitEnergy(): " <<
 274                 fUnitArray[idTPC].GetUnitEnergy() << endl;
 275               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 276               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 277             }
 278
 279             // Fill energy in TPC acceptance
 280             fUnitArray[idTPC].SetUnitEnergy(unitEnergy + ptT[i]);
 281
 282             if(fDebug > 10){
 283               cout << "ptT[" << i << "] : " << ptT[i] << endl;
 284               cout << "unitEnergy in TPC acceptance after : " <<
 285                 fUnitArray[idTPC].GetUnitEnergy() << endl;
 286             }
 287
 288             // Pt cut flag
 289             if(fUnitArray[idTPC].GetUnitEnergy()<ptMin)
 290               fUnitArray[idTPC].SetUnitCutFlag(kPtSmaller);
 291             else fUnitArray[idTPC].SetUnitCutFlag(kPtHigher);
 292             // Detector flag
 293             fUnitArray[idTPC].SetUnitDetectorFlag(kTpc);
 294         }
 295
 296       if(fGrid==1)
 297         {
 298           // Fill track information in EMCAL acceptance
 299           if((etaT[i] >= fEtaMin && etaT[i] <= fEtaMax) &&
 300              (phiT[i] >= fPhiMin && phiT[i] <= fPhiMax))// &&
 301             //   GetCutFlag(i) == 1)
 302             //   ptT[i] > ptMin)
 303             {
 304               trackInEmcalAcc += 1;
 305
 306               if(fDebug>20){
 307                 cout << "before : " << endl;
 308                 cout << "etaT[i] : " << etaT[i] << endl;
 309                 cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
 310               }
 311
 312               // This function should be modified soon
 313 //            Int_t towerID = fGeom->TowerIndexFromEtaPhi(etaT[i],180.0/TMath::Pi()*phiT[i]); // Obsolete
 314               Int_t towerID = fGeom->TowerIndexFromEtaPhi2(etaT[i],180.0/TMath::Pi()*phiT[i]); // Mine modified
 315 //            Int_t towerID = fEMCalGrid->GetIndexFromEtaPhi(phiT[i],etaT[i]);  // Using an EMCal grid -> calculated
 316 //            Int_t towerID = fEMCalGrid->GetIndex(phiT[i],etaT[i]);  // Using an EMCal grid -> tabulated (faster)
 317
 318               Float_t unitEnergy = fUnitArray[towerID].GetUnitEnergy();
 319
 320               if(fDebug>20) {
 321                 cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 322                 cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 323                 cout << "after : " << endl;
 324                 cout << "towerID : " << towerID << endl;
 325                 cout << "etaT[i] : " << etaT[i] << endl;
 326                 cout << "phiT[i](rad) : " << phiT[i] << endl;
 327                 cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
 328                 cout << "unitEnergy in emcal acceptance : " << unitEnergy << endl;
 329                 cout << "fHadCorr : " << fHadCorr << endl;
 330                 cout << "fHCorrection : " << fHCorrection << endl;
 331                 cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 332                 cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 333               }
 334
 335               //OLD WAY:   //Do Hadron Correction
 336               if (fHCorrection != 0)
 337                 {
 338                   Float_t   hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
 339                   unitEnergy -= hCEnergy*TMath::Sin(thetaT[i]);
 340                   if(fDebug>20){
 341                     cout << "Inside loop for hadron correction ==========================" << endl;
 342                     cout << "enT[" << i << "] : " << enT[i] << endl;
 343                     cout << "ptT[" << i << "] : " << ptT[i] << endl;
 344                     cout << "etaT[" << i << "] : " << etaT[i] << endl;
 345                     cout << "Energy correction : " << hCEnergy << endl;
 346                     cout << "unitEnergy : " << unitEnergy << endl;
 347                     cout << "fUnitArray[towerID].GetUnitEnergy() : " <<
 348                       fUnitArray[towerID].GetUnitEnergy() << endl;
 349                   }
 350                 } //end Hadron Correction loop
 351
 352               fUnitArray[towerID].SetUnitEnergy(unitEnergy + ptT[i]);
 353
 354               // Put a pt cut flag
 355               if(fUnitArray[towerID].GetUnitEnergy()<ptMin)
 356                 fUnitArray[towerID].SetUnitCutFlag(kPtSmaller);
 357               else fUnitArray[towerID].SetUnitCutFlag(kPtHigher);
 358               // Detector flag
 359               fUnitArray[towerID].SetUnitDetectorFlag(kTpc);
 360
 361               if(fDebug>10){
 362                 cout << "After pT filled ===============================================" << endl;
 363                 cout << "PtCut : " << ptMin << endl;
 364                 cout << "ptT[" << i << "] : " << ptT[i] << endl;
 365                 cout << "unitEnergy : " << unitEnergy << endl;
 366                 cout << "fUnitArray[towerID].GetUnitEnergy() : " << fUnitArray[towerID].GetUnitEnergy() << endl;
 367                 cout << "fUnitArray[towerID].GetUnitCutFlag() : " << fUnitArray[towerID].GetUnitCutFlag() << endl;
 368                 cout << "fUnitArray[towerID].GetUnitEta() : " << fUnitArray[towerID].GetUnitEta() << endl;
 369                 cout << "fUnitArray[towerID].GetUnitPhi() : " << fUnitArray[towerID].GetUnitPhi() << endl;
 370               }
 371             } // end loop on EMCal acceptance cut + tracks quality
 372           else{
 373             // Outside EMCal acceptance
 374             if(fDebug>2)
 375               cout << "Outside EMCal acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
 376
 377             trackInTpcAcc += 1;
 378
 379             //  Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
 380             Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
 381
 382             Float_t unitEnergy2 = 0.;
 383             unitEnergy2 = fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
 384
 385             if(unitEnergy2 > 0. && fDebug >10){
 386               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 387               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 388               cout << "TPC id : " << idTPC << endl;
 389               cout << "Total id : " << fNumUnits-1+idTPC << endl;
 390               cout << "etaT[" << i << "] : " << etaT[i] << endl;
 391               cout << "phiT[" << i << "] : " << phiT[i] << endl;
 392               cout << "unitEnergy outside emcal acceptance : " << unitEnergy2 << endl;
 393               cout << "fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(): " <<
 394                 fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy() << endl;
 395               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 396               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 397             }
 398
 399             // Fill energy outside emcal acceptance
 400             fUnitArray[fNumUnits-1+idTPC].SetUnitEnergy(unitEnergy2 + ptT[i]);
 401
 402             // Pt cut flag
 403             if(fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy()<ptMin)
 404               fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtSmaller);
 405             else fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtHigher);
 406             // Detector flag
 407             fUnitArray[fNumUnits-1+idTPC].SetUnitDetectorFlag(kTpc);
 408           }
 409         } // end fGrid==1
 410     } // end loop on entries (tpc tracks)
 411
 412   if(fDebug>3){
 413     printf("Number of tracks in EMCAL acceptance: %d\n", trackInEmcalAcc);
 414     printf("Number of tracks in TPC acceptance: %d\n", trackInTpcAcc);
 415   }
 416
 417   if(fGrid==0) fNIn = trackInTpcAccOnly;
 418   if(fGrid==1) fNIn = trackInEmcalAcc+trackInTpcAcc;
 419
 420   fOpt = fReaderHeader->GetDetector();
 421
 422   if(fDebug>1) printf("fOpt in FillUnitArrayFromTPCTracks : %d\n", fOpt);
 423
 424       if(fOpt==1 || option=="tpc") // if only TPC
 425         { //Set all unit flags, Eta, Phi
 426
 427           if(fGrid==0)
 428             {
 429               if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
 430               for(Int_t i=0; i<nElements; i++)
 431                 {
 432                   Float_t eta = 0.;
 433                   Float_t phi = 0.;
 434                   if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
 435                   fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
 436                   fUnitArray[i].SetUnitEntries(nElements);
 437                   fUnitArray[i].SetUnitID(i);
 438                   fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta);
 439                   fUnitArray[i].SetUnitEta(eta);
 440                   fUnitArray[i].SetUnitPhi(phi);
 441                   fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
 442                   fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
 443                 }
 444             }
 445
 446           if(fGrid==1)
 447             {
 448               if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
 449               for(Int_t i=0; i<fNumUnits+nElements; i++)
 450                 {
 451                   Float_t eta = 0.;
 452                   Float_t phi = 0.;
 453                   if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
 454                   //Set all units to be outside a jet initially
 455                   fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
 456                   fUnitArray[i].SetUnitEntries(fNumUnits+nElements);
 457                   fUnitArray[i].SetUnitID(i);
 458
 459                   if(fUnitArray[i].GetUnitID()<fNumUnits)
 460                     {
 461                       // fGeom->EtaPhiFromIndex2(fUnitArray[i].GetUnitID(), eta, phi); // My function in HEADPythia63
 462                       fGeom->EtaPhiFromIndex(fUnitArray[i].GetUnitID(), eta, phi); // From EMCal geometry
 463                       phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 464                       // fEMCalGrid->GetEtaPhiFromIndex2(i,phi,eta); // My function from Grid
 465                       // fEMCalGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta); // My function from Grid
 466                       fUnitArray[i].SetUnitEta(eta);
 467                       //fUnitArray[i].SetUnitPhi(phi*TMath::Pi()/180.0);
 468                       fUnitArray[i].SetUnitPhi(phi);
 469                       fUnitArray[i].SetUnitDeta(fEMCalGrid->GetDeta());
 470                       fUnitArray[i].SetUnitDphi(fEMCalGrid->GetDphi());
 471                     }
 472                   else {
 473                     fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID()+1-fNumUnits,phi,eta);
 474                     fUnitArray[i].SetUnitEta(eta);
 475                     fUnitArray[i].SetUnitPhi(phi);
 476                     fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
 477                     fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
 478
 479                     if(fDebug>10)
 480                       {
 481                         if(fUnitArray[i].GetUnitEnergy()!=0.){
 482                           cout << "(fUnitArray[" << i << "].GetUnitID()+1-fNumUnits : " <<
 483                             fUnitArray[i].GetUnitID()+1-fNumUnits << endl;
 484                           cout << "(fUnitArray[" << i << "].GetUnitEnergy() : " <<
 485                             fUnitArray[i].GetUnitEnergy() << endl;
 486                           cout << "(fUnitArray[" << i << "].GetUnitEta() : " <<
 487                             fUnitArray[i].GetUnitEta() << endl;
 488                           cout << "(fUnitArray[" << i << "].GetUnitPhi() : " <<
 489                             fUnitArray[i].GetUnitPhi() << endl;
 490                         }
 491                       }
 492                   }
 493                   fUnitArray[i].SetUnitClusterID(0);
 494                 }//end loop over all units in array (same as all towers in EMCAL)
 495             }
 496
 497           if(fDebug>20)
 498             {
 499               for(Int_t i=0; i<fNumUnits+nElements; i++)
 500                 {
 501                   if(fUnitArray[i].GetUnitEnergy()!=0.){
 502                     cout << "######################################################### " << endl;
 503                     cout << "Final UnitArray filled with energy != 0" << i << endl;
 504                     cout << "Pointeur UnitArray : " << fUnitArray << " ID : " <<
 505                       fUnitArray[i].GetUnitID() << " Energy : " << fUnitArray[i].GetUnitEnergy() <<
 506                       " Eta : " << fUnitArray[i].GetUnitEta() << " Phi : " << fUnitArray[i].GetUnitPhi() << endl;
 507                   }
 508                 }
 509             }
 510
 511         } // end  if(fOpt==1 || option=="tpc")
 512
 513       delete[] enT;
 514       delete[] ptT;
 515       delete[] etaT;
 516       delete[] phiT;
 517       delete[] thetaT;
 518
 519 }
 520
 521 //__________________________________________________________
 522 void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
 523 {
 524   for(Int_t j=0; j<fNphi+1; j++) {
 525     for(Int_t i=0; i<fNeta+1; i++) {
 526
 527       // TPC grid only
 528       //-------------------------------------
 529       if(fGrid==0) {
 530         if(j*(fNeta+1)+i == index) {
 531           eta = fEta2->At(i);
 532           phi = fPhi2->At(j);
 533         }
 534       }
 535
 536       // TPC-EMCAL grid
 537       //-------------------------------------
 538       Int_t ii = 0;
 539       if(i==0) ii = 0;
 540       if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
 541       if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
 542
 543       if(fGrid==1) {
 544         if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
 545           eta = fEta2->At(i);
 546           phi = fPhi2->At(j);
 547         }
 548
 549         if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
 550            ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
 551           if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
 552           else eta = fEta2->At(i);
 553           phi = fPhi2->At(j);
 554         }
 555
 556         if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
 557           eta = fEta2->At(i);
 558           phi = fPhi2->At(j);
 559         }
 560       }
 561     }
 562   }
 563 }
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575