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