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