Class moved to STEERBase.
[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 "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