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