Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / JETAN / AliJetAODFillUnitArrayTracks.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 //======================================================================
18 // ***July 2006
19 // Fill Unit Array class 
20 // Class used by AliJetESDReader to fill a UnitArray from the information extracted 
21 // from the particle tracks
22 // Author: magali.estienne@ires.in2p3.fr
23 //======================================================================
24
25
26 // --- ROOT system ---
27 #include <TVector3.h>
28 #include <TProcessID.h>
29
30 #include "AliJetUnitArray.h"
31 #include "AliJetHadronCorrectionv1.h"
32 #include "AliJetAODFillUnitArrayTracks.h"
33
34 // --- ROOT system ---
35 class TSystem;
36 class TLorentzVector;
37 class TGeoManager;
38
39 // --- AliRoot header files ---
40 class AliJetFinder;
41
42 using std::cout;
43 using std::endl;
44 ClassImp(AliJetAODFillUnitArrayTracks)
45
46 //_____________________________________________________________________________
47 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks()
48   : AliJetFillUnitArray(),
49     fNumUnits(0),
50     fHadCorr(0),
51     fApplyMIPCorrection(kTRUE),
52     fAOD(0x0),
53     fGrid0(0x0),
54     fGrid1(0x0),
55     fGrid2(0x0),
56     fGrid3(0x0),
57     fGrid4(0x0)
58 {
59   // constructor
60 }
61
62 //_____________________________________________________________________________
63 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(AliAODEvent* aod)
64   : AliJetFillUnitArray(),
65     fNumUnits(0),
66     fHadCorr(0),
67     fApplyMIPCorrection(kTRUE),
68     fAOD(aod),
69     fGrid0(0x0),
70     fGrid1(0x0),
71     fGrid2(0x0),
72     fGrid3(0x0),
73     fGrid4(0x0)
74 {
75   // constructor
76 }
77
78 //_____________________________________________________________________________
79 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(const AliJetAODFillUnitArrayTracks &det)
80   : AliJetFillUnitArray(det),
81     fNumUnits(det.fNumUnits),
82     fHadCorr(det.fHadCorr),
83     fApplyMIPCorrection(det.fApplyMIPCorrection),
84     fAOD(det.fAOD),
85     fGrid0(det.fGrid0),
86     fGrid1(det.fGrid1),
87     fGrid2(det.fGrid2),
88     fGrid3(det.fGrid3),
89     fGrid4(det.fGrid4)
90 {
91   // Copy constructor
92 }
93
94 //_____________________________________________________________________________
95 AliJetAODFillUnitArrayTracks& AliJetAODFillUnitArrayTracks::operator=(const AliJetAODFillUnitArrayTracks& other)
96 {
97   // Assignment
98
99     fNumUnits = other.fNumUnits;
100     fHadCorr = other.fHadCorr;
101     fApplyMIPCorrection = other.fApplyMIPCorrection;
102     fAOD = other.fAOD;
103     fGrid0 = other.fGrid0;
104     fGrid1 = other.fGrid1;
105     fGrid2 = other.fGrid2;
106     fGrid3 = other.fGrid3;
107     fGrid4 = other.fGrid4;
108
109     return (*this);
110 }
111
112 //____________________________________________________________________________
113 void AliJetAODFillUnitArrayTracks::InitParameters()
114 {
115   //  fHadCorr        = 0;     // For hadron correction
116   fNumUnits = fGeom->GetEMCGeometry()->GetNCells();      // Number of towers in EMCAL
117
118   fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, 
119                         fPhiMax,fEtaMin,fEtaMax);
120   fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc, 
121                         fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
122
123   fIndex = fTPCGrid->GetIndexObject();
124
125   if(fDebug>20){
126     for(Int_t i=0; i<fNphi+1; i++)
127       for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
128           (*fIndex)(i,j) << endl; }
129   } 
130   if(fDebug>1) printf("\n Parameters initiated ! \n");
131 }
132
133 //_____________________________________________________________________________
134 AliJetAODFillUnitArrayTracks::~AliJetAODFillUnitArrayTracks()
135 {
136   // destructor
137 }
138
139 //_____________________________________________________________________________
140 void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/)
141 //void AliJetAODFillUnitArrayTracks::Exec(Option_t* option)
142 {
143   //
144   // Main method.
145   // Fill the unit array with the charged particle information in AOD
146   //
147
148   fDebug = fReaderHeader->GetDebug();
149   
150   // Set parameters
151   InitParameters();
152   fRef->Clear();
153
154   // get number of tracks in event (for the loop)
155   Int_t goodTrack = 0;
156   Int_t nt = 0;
157 //(not used ?)  Int_t nt2 = 0;
158   Int_t nmax = 0;
159   Float_t pt,eta,phi;
160   //  Int_t sflag = 0;
161   TVector3 p3;
162
163   nt = fAOD->GetNumberOfTracks();
164   if(fDebug>1) cout << "Number of Tracks in AOD : " << nt << endl;
165
166   // temporary storage of signal and pt cut flag
167   Int_t* sflag  = new Int_t[nt];
168   Int_t* cflag  = new Int_t[nt];
169
170   // get cuts set by user
171   Float_t ptMin  = fReaderHeader->GetPtCut();
172   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
173   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
174   fOpt = fReaderHeader->GetDetector();
175   fDZ  = fReaderHeader->GetDZ();
176   UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
177
178   Int_t nTracksEmcal      = 0;
179   Int_t nTracksEmcalDZ    = 0;
180   Int_t nTracksTpc        = 0;
181   Int_t nTracksTpcOnly    = 0;
182   Int_t nTracksEmcalCut   = 0;
183   Int_t nTracksEmcalDZCut = 0;
184   Int_t nTracksTpcCut     = 0;
185   Int_t nTracksTpcOnlyCut = 0;
186
187   fGrid = fTPCGrid->GetGridType();
188
189   //loop over tracks
190     nmax = nt;  
191     for (Int_t it = 0; it < nmax; it++) {
192       AliAODTrack *track;
193       track = fAOD->GetTrack(it);
194       UInt_t status = track->GetStatus();
195     
196       Double_t mom[3];
197       track->GetPxPyPz(mom);
198
199       p3.SetXYZ(mom[0],mom[1],mom[2]);
200       pt = p3.Pt();
201
202       eta = p3.Eta();
203       phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
204       if (status == 0) continue;
205       if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
206
207       if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
208       
209       if (goodTrack == 0) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
210     
211       sflag[goodTrack]=0;
212       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
213       cflag[goodTrack]=0;
214       if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
215       fRef->Add(track);
216
217       if(fGrid==0)
218         {
219           // Only TPC filled from its grid in its total acceptance
220           
221           Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
222           Bool_t ok = kFALSE;
223           Bool_t ref = kFALSE;
224
225           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
226           TRefArray *reference = uArray->GetUnitTrackRef();
227           if (reference->GetEntries() == 0)  {
228               new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
229           }
230               
231           reference->Add(track);
232
233           Float_t unitEnergy = 0.;
234           unitEnergy = uArray->GetUnitEnergy();
235           // nTracksTpcOnly is necessary to count the number of candidate cells
236           // but it doesn't give the real multiplicity -> it will be extracted 
237           // in the jet finder through the reference to tracks
238           if(unitEnergy==0.){
239             nTracksTpcOnly++;
240             ok = kTRUE;
241             ref = kTRUE;
242           }
243
244           // Fill energy in TPC acceptance
245           uArray->SetUnitEnergy(unitEnergy + pt);
246
247           // Pt cut flag
248           if(uArray->GetUnitEnergy()<ptMin){
249             uArray->SetUnitCutFlag(kPtSmaller);
250           }
251           else {
252             uArray->SetUnitCutFlag(kPtHigher);
253             if(ok) nTracksTpcOnlyCut++;
254           }
255
256           // Signal flag
257           if(sflag[goodTrack] == 1) { 
258             uArray->SetUnitSignalFlag(kGood);
259             uArray->SetUnitSignalFlagC(kFALSE,kGood);
260           } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
261
262           if(uArray->GetUnitEnergy()>0 && ref){
263             if(!fProcId) {
264               fProcId = kTRUE;
265               new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
266             }
267             fRefArray->Add(uArray);
268           }
269         }
270     
271       if(fGrid==1)
272         {
273
274           Int_t nElements = fTPCGrid->GetNEntries();
275           // Fill track information in EMCAL acceptance
276           if((eta >= fEtaMin && eta <= fEtaMax) &&
277              (phi >= fPhiMin && phi <= fPhiMax))// &&
278             {
279               // Include dead-zones
280               if(fDZ)
281                 {
282                   Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
283                   Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
284                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
285                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
286                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
287                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
288                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
289                   Int_t n0 = fGrid0->GetNEntries();
290                   Int_t n1 = fGrid1->GetNEntries();
291                   Int_t n2 = fGrid2->GetNEntries();
292                   Int_t n3 = fGrid3->GetNEntries();
293                   
294                   if(phi >= phimin0 && phi <= phimax0){
295                     Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
296                     AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
297                     TRefArray *reference0 = uArray0->GetUnitTrackRef();
298                     if (reference0->GetEntries() == 0) {
299                         new(reference0) TRefArray(TProcessID::GetProcessWithUID(track));
300                     }
301
302                     reference0->Add(track);
303                     uArray0->SetUnitTrackID(it);
304
305                     Float_t uEnergy0 = uArray0->GetUnitEnergy();
306                     Bool_t ok0 = kFALSE;
307                     Bool_t ref0 = kFALSE;
308                     if(uEnergy0==0.){
309                       nTracksEmcalDZ++;
310                       ok0 = kTRUE;
311                       ref0 = kTRUE;
312                     }
313                     uArray0->SetUnitEnergy(uEnergy0+pt);
314
315                     if(uArray0->GetUnitEnergy()<ptMin)
316                       uArray0->SetUnitCutFlag(kPtSmaller);
317                     else {
318                       uArray0->SetUnitCutFlag(kPtHigher);
319                       if(ok0) nTracksEmcalDZCut++;
320                     }
321                     if(sflag[goodTrack] == 1) {
322                       uArray0->SetUnitSignalFlag(kGood);
323                       uArray0->SetUnitSignalFlagC(kFALSE,kGood);
324                     } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
325
326                     if(uArray0->GetUnitEnergy()>0 && ref0){
327                       if(!fProcId){
328                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
329                         fProcId = kTRUE;
330                       }
331                       fRefArray->Add(uArray0);
332                     }
333                   }
334
335                   if(phi >= phimin1 && phi <= phimax1){
336                     Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
337                     AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
338                     TRefArray *reference1 = uArray1->GetUnitTrackRef();
339                     if (reference1->GetEntries() == 0) {
340                         new(reference1) TRefArray(TProcessID::GetProcessWithUID(track));
341                     }
342
343                     reference1->Add(track);
344
345                     Float_t uEnergy1 = uArray1->GetUnitEnergy();
346                     Bool_t ok1 = kFALSE;
347                     Bool_t ref1 = kFALSE;
348                     if(uEnergy1==0.){
349                       nTracksEmcalDZ++;
350                       ok1 = kTRUE;
351                       ref1 = kTRUE;
352                     }
353                     uArray1->SetUnitEnergy(uEnergy1+pt);
354
355                     if(uArray1->GetUnitEnergy()<ptMin)
356                       uArray1->SetUnitCutFlag(kPtSmaller);
357                     else {
358                       uArray1->SetUnitCutFlag(kPtHigher);
359                       if(ok1) nTracksEmcalDZCut++;
360                     }
361                     if(sflag[goodTrack] == 1) {
362                       uArray1->SetUnitSignalFlag(kGood);
363                       uArray1->SetUnitSignalFlagC(kFALSE,kGood);
364                     } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
365
366                     if(uArray1->GetUnitEnergy()>0 && ref1){
367                       if(!fProcId){
368                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
369                         fProcId = kTRUE;
370                       }
371                       fRefArray->Add(uArray1);
372                     }
373                   }
374
375                   if(phi >= phimin2 && phi <= phimax2){
376                     Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
377                     AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
378                     TRefArray *reference2 = uArray2->GetUnitTrackRef();
379                     if (reference2->GetEntries() == 0) {
380                         new(reference2) TRefArray(TProcessID::GetProcessWithUID(track));
381                     }
382
383                     reference2->Add(track);
384
385                     Float_t uEnergy2 = uArray2->GetUnitEnergy();
386                     Bool_t ok2 = kFALSE;
387                     Bool_t ref2 = kFALSE;
388                     if(uEnergy2==0.){
389                       nTracksEmcalDZ++;
390                       ok2 = kTRUE;
391                       ref2 = kTRUE;
392                     }
393                     uArray2->SetUnitEnergy(uEnergy2+pt);
394
395                     if(uArray2->GetUnitEnergy()<ptMin)
396                       uArray2->SetUnitCutFlag(kPtSmaller);
397                     else {
398                       uArray2->SetUnitCutFlag(kPtHigher);
399                       if(ok2) nTracksEmcalDZCut++;
400                     }
401                     if(sflag[goodTrack] == 1) {
402                       uArray2->SetUnitSignalFlag(kGood);
403                       uArray2->SetUnitSignalFlagC(kFALSE,kGood);
404                     } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
405
406                     if(uArray2->GetUnitEnergy()>0 && ref2){
407                       if(!fProcId){
408                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
409                         fProcId = kTRUE;
410                       }
411                       fRefArray->Add(uArray2);
412                     }
413                   }
414
415                   if(phi >= phimin3 && phi <= phimax3){
416                     Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
417                     AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
418                     TRefArray *reference3 = uArray3->GetUnitTrackRef();
419                     if (reference3->GetEntries() == 0) {
420                         new(reference3) TRefArray(TProcessID::GetProcessWithUID(track));
421                     }
422
423                     reference3->Add(track);
424
425                     Float_t uEnergy3 = uArray3->GetUnitEnergy();
426                     Bool_t ok3 = kFALSE;
427                     Bool_t ref3 = kFALSE;
428                     if(uEnergy3==0.){
429                       nTracksEmcalDZ++;
430                       ok3 = kTRUE;
431                       ref3 = kTRUE;
432                     }
433                     uArray3->SetUnitEnergy(uEnergy3+pt);
434
435                     if(uArray3->GetUnitEnergy()<ptMin)
436                       uArray3->SetUnitCutFlag(kPtSmaller);
437                     else {
438                       uArray3->SetUnitCutFlag(kPtHigher);
439                       if(ok3) nTracksEmcalDZCut++;
440                     }
441                     if(sflag[goodTrack] == 1) {
442                       uArray3->SetUnitSignalFlag(kGood);
443                       uArray3->SetUnitSignalFlagC(kFALSE,kGood);
444                     } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
445
446                     if(uArray3->GetUnitEnergy()>0 && ref3){
447                       if(!fProcId){
448                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
449                         fProcId = kTRUE;
450                       }
451                       fRefArray->Add(uArray3);
452                     }
453                   }
454
455                   if(phi >= phimin4 && phi <= phimax4){
456                     Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
457                     AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
458                     TRefArray *reference4 = uArray4->GetUnitTrackRef();
459                     if (reference4->GetEntries() == 0) {
460                         new(reference4) TRefArray(TProcessID::GetProcessWithUID(track));
461                     }
462
463                     reference4->Add(track);
464
465                     Float_t uEnergy4 = uArray4->GetUnitEnergy();
466                     Bool_t ok4 = kFALSE;
467                     Bool_t ref4 = kFALSE;
468                     if(uEnergy4==0.){
469                       nTracksEmcalDZ++;
470                       ok4 = kTRUE;
471                       ref4 = kTRUE;
472                     }
473                     uArray4->SetUnitEnergy(uEnergy4+pt);
474
475                     if(uArray4->GetUnitEnergy()<ptMin)
476                       uArray4->SetUnitCutFlag(kPtSmaller);
477                     else {
478                       uArray4->SetUnitCutFlag(kPtHigher);
479                       if(ok4) nTracksEmcalDZCut++;
480                     }
481                     if(sflag[goodTrack] == 1) {
482                       uArray4->SetUnitSignalFlag(kGood);
483                       uArray4->SetUnitSignalFlagC(kFALSE,kGood);
484                     } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
485
486                     if(uArray4->GetUnitEnergy()>0 && ref4){
487                       if(!fProcId){
488                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
489                         fProcId = kTRUE;
490                       }
491                       fRefArray->Add(uArray4);
492                     }
493                   }
494                 } // end fDZ
495             
496               Int_t towerID = 0;
497               fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
498               if(towerID==-1) continue;
499               
500               AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
501               TRefArray *reference = uArray->GetUnitTrackRef();
502               if (reference->GetEntries() == 0) {
503                   new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
504               }
505
506               reference->Add(track);
507
508               Float_t unitEnergy = uArray->GetUnitEnergy(); 
509               Bool_t ok = kFALSE;
510               Bool_t ref = kFALSE;
511               if(unitEnergy==0.){
512                 nTracksEmcal++;
513                 ok=kTRUE;
514                 ref=kTRUE;
515               }
516
517               // Do Hadron Correction
518               // For the moment I apply MIP correction if p >= 0.5 GeV/c
519               // and only for tracks inside EMCal acceptance
520               // End of if(fGrid==1) -> hadron correction for all tracks
521 // Temporarily commented - will be updated removing the AliAODpid dependence
522 /*
523               if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5) 
524                 {
525                   ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETEV1",1.);
526
527                   Double_t etaOut = 0.;
528                   Double_t phiOut = 0.;
529                   ((AliJetHadronCorrectionv1*)fHadCorr)->TrackPositionEMCal(track,etaOut,phiOut);
530
531                   // One has to make sure that the track hits the calorimeter
532                   // We can then correct on average for these particles
533                   if((etaOut >= fEtaMin && etaOut <= fEtaMax) &&
534                      (phiOut >= fPhiMin && phiOut <= fPhiMax))// &&
535                     {
536                       Double_t   hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
537                       unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
538                     } 
539                 } //end Hadron Correction loop
540 */
541
542               if((unitEnergy + pt) > 0.) uArray->SetUnitEnergy(unitEnergy + pt);
543               else uArray->SetUnitEnergy(0.);
544
545               // Put a pt cut flag
546               if(uArray->GetUnitEnergy()<ptMin){
547                 uArray->SetUnitCutFlag(kPtSmaller);
548               }
549               else {
550                 uArray->SetUnitCutFlag(kPtHigher);
551                 if(ok) nTracksEmcalCut++;
552               }
553
554               // Signal flag
555               if(sflag[goodTrack] == 1) { 
556                 uArray->SetUnitSignalFlag(kGood);
557                 uArray->SetUnitSignalFlagC(kFALSE,kGood);
558               } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
559
560               if(uArray->GetUnitEnergy()>0 && ref){
561                 if(!fProcId){
562                   new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
563                   fProcId = kTRUE;
564                 }
565                 fRefArray->Add(uArray);
566               }
567             } // end loop on EMCal acceptance cut + tracks quality
568           else{ 
569             // Outside EMCal acceptance
570             //  Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
571             Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
572
573             AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
574             TRefArray *reference = uArray->GetUnitTrackRef();
575             if (reference->GetEntries() == 0) {
576                 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
577             }
578             reference->Add(track);
579
580             Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
581             Bool_t okout = kFALSE;
582             Bool_t refout = kFALSE;
583             if(unitEnergy2==0.){
584               nTracksTpc++;
585               okout=kTRUE;
586               refout=kTRUE;
587             }
588             // Fill energy outside emcal acceptance
589             uArray->SetUnitEnergy(unitEnergy2 + pt);
590           
591             // Pt cut flag
592             if(uArray->GetUnitEnergy()<ptMin){
593               uArray->SetUnitCutFlag(kPtSmaller);
594             }
595             else {
596               uArray->SetUnitCutFlag(kPtHigher);
597               if(okout) nTracksTpcCut++;
598             }
599             // Signal flag
600             if(sflag[goodTrack] == 1) {
601               uArray->SetUnitSignalFlag(kGood);
602               uArray->SetUnitSignalFlagC(kFALSE,kGood);
603             } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
604
605             if(uArray->GetUnitEnergy()>0 && refout){
606               if(!fProcId){
607                 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
608                 fProcId = kTRUE;
609               }
610               fRefArray->Add(uArray);
611             }
612           }
613
614           /*
615           // Do hadron correction
616           // In principle, the best thing to do as particles which
617           // eta,phi (at vertex) correspond to EMCal dead zone can deposit energy
618           // in the calorimeter as well as particles with eta,phi (at vertex)
619           // outside the EMCal acceptance
620           cout << "Hadronic correction for all tracks goes here" << endl;
621           
622           // For the moment I apply MIP correction if p >= 0.5 GeV/c
623           if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5) 
624             { 
625               ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETEV1",1.);
626               
627               Double_t etaOut = 0.;
628               Double_t phiOut = 0.;
629               Int_t towerID2 = 0;
630               Float_t unitEnergy = 0.;
631               ((AliJetHadronCorrectionv1*)fHadCorr)->TrackPositionEMCal(track,etaOut,phiOut);
632               
633               // One has to make sure that the track hits the calorimeter
634               // We can then correct on average for these particles
635               if((etaOut >= fEtaMin && etaOut <= fEtaMax) &&
636                  (phiOut >= fPhiMin && phiOut <= fPhiMax))// &&
637                 {
638                   Int_t towerID = 0;
639                   Int_t towerID2 = 0;
640                   fGeom->GetAbsCellIdFromEtaPhi(etaOut,phiOut,towerID2);
641                   if(towerID2==-1) continue;
642                   
643                   AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID2);
644                   unitEnergy = uArray->GetUnitEnergy(); 
645                   Bool_t ok = kFALSE;
646                   Bool_t ref = kFALSE;
647                   if(unitEnergy==0.){
648                     nTracksEmcal++;
649                     ok=kTRUE;
650                     ref=kTRUE;
651                   }
652                   
653                   Double_t   hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
654                   unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
655                   
656                   // The energy in this unitarray can be negative
657                   // If the Digits are then filled, make sure the energy sum
658                   // is positive. If not, put the value to zero in AliJetAODUnitArrayEMCalDigit.
659                   // If they are not filled, put this value to zero
660                   uArray->SetUnitEnergy(unitEnergy);
661                   
662                   if(uArray->GetUnitEnergy()!=0. && ref){
663                     if(!fProcId){
664                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
665                       fProcId = kTRUE;
666                     }
667                     fRefArray->Add(uArray);
668                   }
669                 }
670               else cout << "Track out of EMCal acceptance" << endl; 
671               
672             } //end Hadron Correction loop
673           */
674           
675         } // end fGrid==1
676
677       goodTrack++;
678
679     } // end loop on entries (tpc tracks)
680
681
682   if(fDebug>0)
683     {
684       cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
685       cout << "goodTracks: " << goodTrack << endl;
686     }
687
688   fSignalFlag.Set(goodTrack,sflag);
689   fCutFlag.Set(goodTrack,cflag);
690
691   delete[] sflag;
692   delete[] cflag;
693
694   //    } // end loop on entries (tpc tracks)
695       
696   if(fGrid==0) {
697     fNTracks = nTracksTpcOnly;
698     fNTracksCut = nTracksTpcOnlyCut;
699     if(fDebug>10){
700       cout << "fNTracks : " << fNTracks << endl;
701       cout << "fNTracksCut : " << fNTracksCut << endl;
702     }
703   }
704   if(fGrid==1) {
705     fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
706     fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
707     if(fDebug>10){
708       cout << "fNTracks : " << fNTracks << endl;
709       cout << "fNTracksCut : " << fNTracksCut << endl;
710     }
711   }  
712   
713 }
714