JETAN/AliJetESDFillUnitArrayTracks.cxx

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