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