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