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