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 "AliJetGrid.h"
  74
  75 ClassImp(AliJetFillUnitArrayTracks)
  76
  77 //_____________________________________________________________________________
  78 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
  79     : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
  80       fNumUnits(0),
  81       fEtaMinCal(0.),
  82       fEtaMaxCal(0.),
  83       fPhiMinCal(0.),
  84       fPhiMaxCal(0.),
  85       fHadCorr(0x0),
  86       fHCorrection(0),
  87       fNIn(0),
  88       fOpt(0),
  89       fDebug(0),
  90       fReaderHeader(0x0),
  91       fMomentumArray(0x0),
  92       fUnitArray(0x0),
  93       fTPCGrid(0x0),
  94       fEMCalGrid(0x0),
  95       fGeom(0x0),
  96       fNphi(0),
  97       fNeta(0),
  98       fPhi2(0x0),
  99       fEta2(0x0),
 100       fPhi(0x0),
 101       fEta(0x0),
 102       fIndex(0x0),
 103       fParams(0x0),
 104       fGrid(0),
 105       fPhiMin(0.),
 106       fPhiMax(0.),
 107       fEtaMin(0.),
 108       fEtaMax(0.),
 109       fEtaBinInTPCAcc(0),
 110       fPhiBinInTPCAcc(0),
 111       fEtaBinInEMCalAcc(0),
 112       fPhiBinInEMCalAcc(0),
 113       fNbinPhi(0)
 114 {
 115   // constructor
 116 }
 117
 118 //____________________________________________________________________________
 119 void AliJetFillUnitArrayTracks::SetEMCALGeometry()
 120 {
 121     // Set EMCAL geometry information
 122     fGeom = AliJetDummyGeo::GetInstance();
 123     if (fGeom == 0)
 124         fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
 125     if(fDebug>1) printf("\n EMCAL Geometry setted ! \n");
 126 }
 127
 128 //____________________________________________________________________________
 129 void AliJetFillUnitArrayTracks::InitParameters()
 130 {
 131   fHCorrection    = 0;     // For hadron correction
 132   fHadCorr        = 0;     // For hadron correction
 133   fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
 134   fDebug = fReaderHeader->GetDebug();
 135
 136   fEtaMinCal = fGeom->GetArm1EtaMin();
 137   fEtaMaxCal = fGeom->GetArm1EtaMax();
 138   fPhiMinCal = fGeom->GetArm1PhiMin();
 139   fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
 140
 141   if(fDebug>30){
 142     cout << "fEtaMinCal : " << fEtaMinCal << endl;
 143     cout << "fEtaMaxCal : " << fEtaMaxCal << endl;
 144     cout << "fPhiMinCal : " << fPhiMinCal << endl;
 145     cout << "fPhiMaxCal : " << fPhiMaxCal << endl;
 146   }
 147
 148   fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
 149                         fPhiMax,fEtaMin,fEtaMax);
 150   fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
 151                         fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
 152
 153   fEta   = fTPCGrid->GetArrayEta();
 154   fPhi   = fTPCGrid->GetArrayPhi();
 155   fIndex = fTPCGrid->GetIndexObject();
 156
 157   if(fDebug>20){
 158     for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
 159     for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
 160
 161     for(Int_t i=0; i<fNphi+1; i++)
 162       for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
 163           (*fIndex)(i,j) << endl; }
 164   }
 165   if(fDebug>1) printf("\n Parameters initiated ! \n");
 166 }
 167
 168 //_____________________________________________________________________________
 169 AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
 170 {
 171   // destructor
 172 }
 173
 174 //_____________________________________________________________________________
 175 void AliJetFillUnitArrayTracks::Exec(Option_t* option)
 176 {
 177   //
 178   // Main method.
 179   // Explain
 180
 181   fDebug = fReaderHeader->GetDebug();
 182   if(fDebug>1) printf("In AliJetFillUnitArrayTracks !");
 183   if(fDebug>3) printf("\nfGeom->GetEntries() = %d\n", fGeom->GetNCells());
 184   // Set EMCal Geometry
 185   SetEMCALGeometry();
 186   // Set parameters
 187   InitParameters();
 188
 189   TClonesArray *lvArray = fMomentumArray;     // Correct checked !
 190   Int_t nInT =  lvArray->GetEntries();        // Correct checked !
 191   Float_t ptMin = fReaderHeader->GetPtCut();  // Correct checked !
 192
 193   // sflag -> not yet implemented !!!!
 194
 195   if(fDebug>3) cout << "nInT : " << nInT << endl;
 196
 197   if (nInT == 0) return;
 198
 199   // local arrays for input
 200   Float_t* enT  = new Float_t[nInT];
 201   Float_t* ptT  = new Float_t[nInT];
 202   Float_t* etaT = new Float_t[nInT];
 203   Float_t* phiT = new Float_t[nInT];
 204   Float_t* thetaT = new Float_t[nInT];
 205
 206   Int_t trackInEmcalAcc = 0;
 207   Int_t trackInTpcAcc = 0;
 208   Int_t trackInTpcAccOnly = 0;
 209
 210   Int_t nElements = fTPCGrid->GetNEntries();
 211   Int_t nElements2 = fEMCalGrid->GetNEntries();
 212   fGrid = fTPCGrid->GetGridType();
 213
 214   if(fDebug>3){
 215     cout << "nElements : " << nElements << endl;
 216     cout << "nElements2 : " << nElements2 << endl;
 217     cout << "fNumUnits : " << fNumUnits << endl;
 218     cout << "sum : " << fNumUnits+nElements << endl;
 219   }
 220
 221   // Set energy exactly to zero
 222   if(fGrid==0)
 223     for(Int_t k=0; k<nElements; k++)
 224       fUnitArray[k].SetUnitEnergy(0.);
 225
 226   if(fGrid==1)
 227     for(Int_t k=0; k<fNumUnits+nElements; k++)
 228       fUnitArray[k].SetUnitEnergy(0.);
 229
 230   // load input vectors
 231   for (Int_t i = 0; i < nInT; i++)
 232     {
 233       TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
 234       enT[i]  = lv->Energy();
 235       ptT[i]  = lv->Pt();
 236       etaT[i] = lv->Eta();
 237       phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
 238       thetaT[i] = 2.0*TMath::ATan(TMath::Exp(-etaT[i]));
 239
 240       if(fDebug>20){
 241         cout << "enT[" << i << "] : " <<  enT[i] << endl;
 242         cout << "ptT[" << i << "] : " <<  ptT[i] << endl;
 243         cout << "etaT[" << i << "] : " <<  etaT[i] << endl;
 244         cout << "phiT[" << i << "] : " <<  phiT[i] << endl;
 245         cout << "thetaT[" << i << "] : " <<  thetaT[i] << endl;
 246         cout << "fEtaMinCal : " << fEtaMinCal << ", fEtaMaxCal : " << fEtaMaxCal << endl;
 247         cout << "fPhiMinCal : " << fPhiMinCal << ", fPhiMaxCal : " << fPhiMaxCal << endl;
 248         cout << "fEtaMin : " << fEtaMin << ", fEtaMax : " << fEtaMax << endl;
 249         cout << "fPhiMin : " << fPhiMin << ", fPhiMax : " << fPhiMax << endl;
 250       }
 251
 252       if(fGrid==0)
 253         {
 254           // For the moment, only TPC filled from its grid in its total acceptance
 255           if(fDebug>2)
 256             cout << "In total TPC acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
 257
 258             trackInTpcAccOnly += 1;
 259
 260             Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
 261
 262             Float_t unitEnergy = 0.;
 263             unitEnergy = fUnitArray[idTPC].GetUnitEnergy();
 264
 265             if(unitEnergy > 0. && fDebug >10){
 266               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 267               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 268               cout << "TPC id : " << idTPC << endl;
 269               cout << "etaT[" << i << "] : " << etaT[i] << endl;
 270               cout << "phiT[" << i << "] : " << phiT[i] << endl;
 271               cout << "unitEnergy in TPC acceptance : " << unitEnergy << endl;
 272               cout << "fUnitArray[idTPC].GetUnitEnergy(): " <<
 273                 fUnitArray[idTPC].GetUnitEnergy() << endl;
 274               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 275               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 276             }
 277
 278             // Fill energy in TPC acceptance
 279             fUnitArray[idTPC].SetUnitEnergy(unitEnergy + ptT[i]);
 280
 281             if(fDebug > 10){
 282               cout << "ptT[" << i << "] : " << ptT[i] << endl;
 283               cout << "unitEnergy in TPC acceptance after : " <<
 284                 fUnitArray[idTPC].GetUnitEnergy() << endl;
 285             }
 286
 287             // Pt cut flag
 288             if(fUnitArray[idTPC].GetUnitEnergy()<ptMin)
 289               fUnitArray[idTPC].SetUnitCutFlag(kPtSmaller);
 290             else fUnitArray[idTPC].SetUnitCutFlag(kPtHigher);
 291             // Detector flag
 292             fUnitArray[idTPC].SetUnitDetectorFlag(kTpc);
 293         }
 294
 295       if(fGrid==1)
 296         {
 297           // Fill track information in EMCAL acceptance
 298           if((etaT[i] >= fEtaMin && etaT[i] <= fEtaMax) &&
 299              (phiT[i] >= fPhiMin && phiT[i] <= fPhiMax))// &&
 300             //   GetCutFlag(i) == 1)
 301             //   ptT[i] > ptMin)
 302             {
 303               trackInEmcalAcc += 1;
 304
 305               if(fDebug>20){
 306                 cout << "before : " << endl;
 307                 cout << "etaT[i] : " << etaT[i] << endl;
 308                 cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
 309               }
 310
 311               // This function should be modified soon
 312 //            Int_t towerID = fGeom->TowerIndexFromEtaPhi(etaT[i],180.0/TMath::Pi()*phiT[i]); // Obsolete
 313               Int_t towerID = fGeom->TowerIndexFromEtaPhi2(etaT[i],180.0/TMath::Pi()*phiT[i]); // Mine modified
 314 //            Int_t towerID = fEMCalGrid->GetIndexFromEtaPhi(phiT[i],etaT[i]);  // Using an EMCal grid -> calculated
 315 //            Int_t towerID = fEMCalGrid->GetIndex(phiT[i],etaT[i]);  // Using an EMCal grid -> tabulated (faster)
 316
 317               Float_t unitEnergy = fUnitArray[towerID].GetUnitEnergy();
 318
 319               if(fDebug>20) {
 320                 cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 321                 cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 322                 cout << "after : " << endl;
 323                 cout << "towerID : " << towerID << endl;
 324                 cout << "etaT[i] : " << etaT[i] << endl;
 325                 cout << "phiT[i](rad) : " << phiT[i] << endl;
 326                 cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
 327                 cout << "unitEnergy in emcal acceptance : " << unitEnergy << endl;
 328                 cout << "fHadCorr : " << fHadCorr << endl;
 329                 cout << "fHCorrection : " << fHCorrection << endl;
 330                 cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 331                 cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 332               }
 333
 334               //OLD WAY:   //Do Hadron Correction
 335               if (fHCorrection != 0)
 336                 {
 337                   Float_t   hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
 338                   unitEnergy -= hCEnergy*TMath::Sin(thetaT[i]);
 339                   if(fDebug>20){
 340                     cout << "Inside loop for hadron correction ==========================" << endl;
 341                     cout << "enT[" << i << "] : " << enT[i] << endl;
 342                     cout << "ptT[" << i << "] : " << ptT[i] << endl;
 343                     cout << "etaT[" << i << "] : " << etaT[i] << endl;
 344                     cout << "Energy correction : " << hCEnergy << endl;
 345                     cout << "unitEnergy : " << unitEnergy << endl;
 346                     cout << "fUnitArray[towerID].GetUnitEnergy() : " <<
 347                       fUnitArray[towerID].GetUnitEnergy() << endl;
 348                   }
 349                 } //end Hadron Correction loop
 350
 351               fUnitArray[towerID].SetUnitEnergy(unitEnergy + ptT[i]);
 352
 353               // Put a pt cut flag
 354               if(fUnitArray[towerID].GetUnitEnergy()<ptMin)
 355                 fUnitArray[towerID].SetUnitCutFlag(kPtSmaller);
 356               else fUnitArray[towerID].SetUnitCutFlag(kPtHigher);
 357               // Detector flag
 358               fUnitArray[towerID].SetUnitDetectorFlag(kTpc);
 359
 360               if(fDebug>10){
 361                 cout << "After pT filled ===============================================" << endl;
 362                 cout << "PtCut : " << ptMin << endl;
 363                 cout << "ptT[" << i << "] : " << ptT[i] << endl;
 364                 cout << "unitEnergy : " << unitEnergy << endl;
 365                 cout << "fUnitArray[towerID].GetUnitEnergy() : " << fUnitArray[towerID].GetUnitEnergy() << endl;
 366                 cout << "fUnitArray[towerID].GetUnitCutFlag() : " << fUnitArray[towerID].GetUnitCutFlag() << endl;
 367                 cout << "fUnitArray[towerID].GetUnitEta() : " << fUnitArray[towerID].GetUnitEta() << endl;
 368                 cout << "fUnitArray[towerID].GetUnitPhi() : " << fUnitArray[towerID].GetUnitPhi() << endl;
 369               }
 370             } // end loop on EMCal acceptance cut + tracks quality
 371           else{
 372             // Outside EMCal acceptance
 373             if(fDebug>2)
 374               cout << "Outside EMCal acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
 375
 376             trackInTpcAcc += 1;
 377
 378             //  Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
 379             Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
 380
 381             Float_t unitEnergy2 = 0.;
 382             unitEnergy2 = fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
 383
 384             if(unitEnergy2 > 0. && fDebug >10){
 385               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 386               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 387               cout << "TPC id : " << idTPC << endl;
 388               cout << "Total id : " << fNumUnits-1+idTPC << endl;
 389               cout << "etaT[" << i << "] : " << etaT[i] << endl;
 390               cout << "phiT[" << i << "] : " << phiT[i] << endl;
 391               cout << "unitEnergy outside emcal acceptance : " << unitEnergy2 << endl;
 392               cout << "fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(): " <<
 393                 fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy() << endl;
 394               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 395               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
 396             }
 397
 398             // Fill energy outside emcal acceptance
 399             fUnitArray[fNumUnits-1+idTPC].SetUnitEnergy(unitEnergy2 + ptT[i]);
 400
 401             // Pt cut flag
 402             if(fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy()<ptMin)
 403               fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtSmaller);
 404             else fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtHigher);
 405             // Detector flag
 406             fUnitArray[fNumUnits-1+idTPC].SetUnitDetectorFlag(kTpc);
 407           }
 408         } // end fGrid==1
 409     } // end loop on entries (tpc tracks)
 410
 411   if(fDebug>3){
 412     printf("Number of tracks in EMCAL acceptance: %d\n", trackInEmcalAcc);
 413     printf("Number of tracks in TPC acceptance: %d\n", trackInTpcAcc);
 414   }
 415
 416   if(fGrid==0) fNIn = trackInTpcAccOnly;
 417   if(fGrid==1) fNIn = trackInEmcalAcc+trackInTpcAcc;
 418
 419   fOpt = fReaderHeader->GetDetector();
 420
 421   if(fDebug>1) printf("fOpt in FillUnitArrayFromTPCTracks : %d\n", fOpt);
 422
 423       if(fOpt==1 || option=="tpc") // if only TPC
 424         { //Set all unit flags, Eta, Phi
 425
 426           if(fGrid==0)
 427             {
 428               if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
 429               for(Int_t i=0; i<nElements; i++)
 430                 {
 431                   Float_t eta = 0.;
 432                   Float_t phi = 0.;
 433                   if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
 434                   fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
 435                   fUnitArray[i].SetUnitEntries(nElements);
 436                   fUnitArray[i].SetUnitID(i);
 437                   fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta);
 438                   fUnitArray[i].SetUnitEta(eta);
 439                   fUnitArray[i].SetUnitPhi(phi);
 440                   fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
 441                   fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
 442                 }
 443             }
 444
 445           if(fGrid==1)
 446             {
 447               if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
 448               for(Int_t i=0; i<fNumUnits+nElements; i++)
 449                 {
 450                   Float_t eta = 0.;
 451                   Float_t phi = 0.;
 452                   if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
 453                   //Set all units to be outside a jet initially
 454                   fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
 455                   fUnitArray[i].SetUnitEntries(fNumUnits+nElements);
 456                   fUnitArray[i].SetUnitID(i);
 457
 458                   if(fUnitArray[i].GetUnitID()<fNumUnits)
 459                     {
 460                       // fGeom->EtaPhiFromIndex2(fUnitArray[i].GetUnitID(), eta, phi); // My function in HEADPythia63
 461                       fGeom->EtaPhiFromIndex(fUnitArray[i].GetUnitID(), eta, phi); // From EMCal geometry
 462                       phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 463                       // fEMCalGrid->GetEtaPhiFromIndex2(i,phi,eta); // My function from Grid
 464                       // fEMCalGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta); // My function from Grid
 465                       fUnitArray[i].SetUnitEta(eta);
 466                       //fUnitArray[i].SetUnitPhi(phi*TMath::Pi()/180.0);
 467                       fUnitArray[i].SetUnitPhi(phi);
 468                       fUnitArray[i].SetUnitDeta(fEMCalGrid->GetDeta());
 469                       fUnitArray[i].SetUnitDphi(fEMCalGrid->GetDphi());
 470                     }
 471                   else {
 472                     fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID()+1-fNumUnits,phi,eta);
 473                     fUnitArray[i].SetUnitEta(eta);
 474                     fUnitArray[i].SetUnitPhi(phi);
 475                     fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
 476                     fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
 477
 478                     if(fDebug>10)
 479                       {
 480                         if(fUnitArray[i].GetUnitEnergy()!=0.){
 481                           cout << "(fUnitArray[" << i << "].GetUnitID()+1-fNumUnits : " <<
 482                             fUnitArray[i].GetUnitID()+1-fNumUnits << endl;
 483                           cout << "(fUnitArray[" << i << "].GetUnitEnergy() : " <<
 484                             fUnitArray[i].GetUnitEnergy() << endl;
 485                           cout << "(fUnitArray[" << i << "].GetUnitEta() : " <<
 486                             fUnitArray[i].GetUnitEta() << endl;
 487                           cout << "(fUnitArray[" << i << "].GetUnitPhi() : " <<
 488                             fUnitArray[i].GetUnitPhi() << endl;
 489                         }
 490                       }
 491                   }
 492                   fUnitArray[i].SetUnitClusterID(0);
 493                 }//end loop over all units in array (same as all towers in EMCAL)
 494             }
 495
 496           if(fDebug>20)
 497             {
 498               for(Int_t i=0; i<fNumUnits+nElements; i++)
 499                 {
 500                   if(fUnitArray[i].GetUnitEnergy()!=0.){
 501                     cout << "######################################################### " << endl;
 502                     cout << "Final UnitArray filled with energy != 0" << i << endl;
 503                     cout << "Pointeur UnitArray : " << fUnitArray << " ID : " <<
 504                       fUnitArray[i].GetUnitID() << " Energy : " << fUnitArray[i].GetUnitEnergy() <<
 505                       " Eta : " << fUnitArray[i].GetUnitEta() << " Phi : " << fUnitArray[i].GetUnitPhi() << endl;
 506                   }
 507                 }
 508             }
 509
 510         } // end  if(fOpt==1 || option=="tpc")
 511
 512       delete[] enT;
 513       delete[] ptT;
 514       delete[] etaT;
 515       delete[] phiT;
 516       delete[] thetaT;
 517
 518 }
 519
 520 //__________________________________________________________
 521 void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
 522 {
 523   for(Int_t j=0; j<fNphi+1; j++) {
 524     for(Int_t i=0; i<fNeta+1; i++) {
 525
 526       // TPC grid only
 527       //-------------------------------------
 528       if(fGrid==0) {
 529         if(j*(fNeta+1)+i == index) {
 530           eta = fEta2->At(i);
 531           phi = fPhi2->At(j);
 532         }
 533       }
 534
 535       // TPC-EMCAL grid
 536       //-------------------------------------
 537       Int_t ii = 0;
 538       if(i==0) ii = 0;
 539       if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
 540       if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
 541
 542       if(fGrid==1) {
 543         if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
 544           eta = fEta2->At(i);
 545           phi = fPhi2->At(j);
 546         }
 547
 548         if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
 549            ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
 550           if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
 551           else eta = fEta2->At(i);
 552           phi = fPhi2->At(j);
 553         }
 554
 555         if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
 556           eta = fEta2->At(i);
 557           phi = fPhi2->At(j);
 558         }
 559       }
 560     }
 561   }
 562 }
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574