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