JETAN/AliJetFillUnitArrayTracks.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 // --- Standard library ---
  27 #include <Riostream.h>
  28
  29 // --- ROOT system ---
  30 #include <TSystem.h>
  31 #include <TLorentzVector.h>
  32 #include <TRefArray.h>
  33 #include <TVector3.h>
  34 #include "TTask.h"
  35 #include <TGeoManager.h>
  36 #include <TMatrixD.h>
  37 #include <TArrayD.h>
  38 #include <TMath.h>
  39 #include <TClonesArray.h>
  40 #include <TProcessID.h>
  41
  42 // --- AliRoot header files ---
  43 #include "AliJetFinder.h"
  44 #include "AliJetReaderHeader.h"
  45 #include "AliJetReader.h"
  46 #include "AliJetESDReader.h"
  47 #include "AliJetESDReaderHeader.h"
  48 #include "AliESDEvent.h"
  49 #include "AliJetDummyGeo.h"
  50 #include "AliJetUnitArray.h"
  51 #include "AliJetFillUnitArrayTracks.h"
  52 #include "AliJetHadronCorrectionv1.h"
  53 #include "AliJetGrid.h"
  54
  55 ClassImp(AliJetFillUnitArrayTracks)
  56
  57 //_____________________________________________________________________________
  58 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
  59   : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
  60     fNumUnits(0),
  61     fEtaMinCal(0),
  62     fEtaMaxCal(0),
  63     fPhiMinCal(0),
  64     fPhiMaxCal(0),
  65     fHadCorr(0),
  66     fHCorrection(0),
  67     fNTracks(0),
  68     fNTracksCut(0),
  69     fOpt(0),
  70     fDZ(0),
  71     fDebug(0),
  72     fReaderHeader(0x0),
  73     fMomentumArray(0x0),
  74     fUnitArray(0x0),
  75     fRefArray(0x0),
  76     fProcId(kFALSE),
  77     fTPCGrid(0x0),
  78     fEMCalGrid(0x0),
  79     fGeom(0x0),
  80     fESD(0x0),
  81     fGrid0(0x0),
  82     fGrid1(0x0),
  83     fGrid2(0x0),
  84     fGrid3(0x0),
  85     fGrid4(0x0),
  86     fNphi(0),
  87     fNeta(0),
  88     fPhi2(0x0),
  89     fEta2(0x0),
  90     fPhi(0x0),
  91     fEta(0x0),
  92     fIndex(0x0),
  93     fParams(0x0),
  94     fGrid(0),
  95     fPhiMin(0),
  96     fPhiMax(0),
  97     fEtaMin(0),
  98     fEtaMax(0),
  99     fEtaBinInTPCAcc(0),
 100     fPhiBinInTPCAcc(0),
 101     fEtaBinInEMCalAcc(0),
 102     fPhiBinInEMCalAcc(0),
 103     fNbinPhi(0)
 104 {
 105   // constructor
 106 }
 107
 108 //_____________________________________________________________________________
 109 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* esd)
 110   : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
 111     fNumUnits(0),
 112     fEtaMinCal(0),
 113     fEtaMaxCal(0),
 114     fPhiMinCal(0),
 115     fPhiMaxCal(0),
 116     fHadCorr(0),
 117     fHCorrection(0),
 118     fNTracks(0),
 119     fNTracksCut(0),
 120     fOpt(0),
 121     fDZ(0),
 122     fDebug(0),
 123     fReaderHeader(0x0),
 124     fMomentumArray(0x0),
 125     fUnitArray(0x0),
 126     fRefArray(0x0),
 127     fProcId(kFALSE),
 128     fTPCGrid(0x0),
 129     fEMCalGrid(0x0),
 130     fGeom(0x0),
 131     fESD(esd),
 132     fGrid0(0x0),
 133     fGrid1(0x0),
 134     fGrid2(0x0),
 135     fGrid3(0x0),
 136     fGrid4(0x0),// first: loop over tracks in ESD
 137     fNphi(0),
 138     fNeta(0),
 139     fPhi2(0x0),
 140     fEta2(0x0),
 141     fPhi(0x0),
 142     fEta(0x0),
 143     fIndex(0x0),
 144     fParams(0x0),
 145     fGrid(0),
 146     fPhiMin(0),
 147     fPhiMax(0),
 148     fEtaMin(0),
 149     fEtaMax(0),
 150     fEtaBinInTPCAcc(0),
 151     fPhiBinInTPCAcc(0),
 152     fEtaBinInEMCalAcc(0),
 153     fPhiBinInEMCalAcc(0),
 154     fNbinPhi(0)
 155 {
 156   // constructor 2
 157 }
 158
 159 //_____________________________________________________________________________
 160 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(const AliJetFillUnitArrayTracks &det)
 161   : TTask(det),//"AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
 162     fNumUnits(det.fNumUnits),
 163     fEtaMinCal(det.fEtaMinCal),
 164     fEtaMaxCal(det.fEtaMaxCal),
 165     fPhiMinCal(det.fPhiMinCal),
 166     fPhiMaxCal(det.fPhiMaxCal),
 167     fHadCorr(det.fHadCorr),
 168     fHCorrection(det.fHCorrection),
 169     fNTracks(det.fNTracks),
 170     fNTracksCut(det.fNTracksCut),
 171     fOpt(det.fOpt),
 172     fDZ(det.fDZ),
 173     fDebug(det.fDebug),
 174     fReaderHeader(det.fReaderHeader),
 175     fMomentumArray(det.fMomentumArray),
 176     fUnitArray(det.fUnitArray),
 177     fRefArray(det.fRefArray),
 178     fProcId(det.fProcId),
 179     fTPCGrid(det.fTPCGrid),
 180     fEMCalGrid(det.fEMCalGrid),
 181     fGeom(det.fGeom),
 182     fESD(det.fESD),
 183     fGrid0(det.fGrid0),
 184     fGrid1(det.fGrid1),
 185     fGrid2(det.fGrid2),
 186     fGrid3(det.fGrid3),
 187     fGrid4(det.fGrid4),
 188     fNphi(det.fNphi),
 189     fNeta(det.fNeta),
 190     fPhi2(det.fPhi2),
 191     fEta2(det.fEta2),
 192     fPhi(det.fPhi),
 193     fEta(det.fEta),
 194     fIndex(det.fIndex),
 195     fParams(det.fParams),
 196     fGrid(det.fGrid),
 197     fPhiMin(det.fPhiMin),
 198     fPhiMax(det.fPhiMax),
 199     fEtaMin(det.fEtaMin),
 200     fEtaMax(det.fEtaMax),
 201     fEtaBinInTPCAcc(det.fEtaBinInTPCAcc),
 202     fPhiBinInTPCAcc(det.fPhiBinInTPCAcc),
 203     fEtaBinInEMCalAcc(det.fEtaBinInEMCalAcc),
 204     fPhiBinInEMCalAcc(det.fPhiBinInEMCalAcc),
 205     fNbinPhi(det.fNbinPhi)
 206 {
 207   // Copy constructor
 208 }
 209
 210 //------------------------------------------------------------------
 211 AliJetFillUnitArrayTracks& AliJetFillUnitArrayTracks::operator=(const AliJetFillUnitArrayTracks& other)
 212 {
 213   // Assignment
 214
 215     fNumUnits = other.fNumUnits;
 216     fEtaMinCal = other.fEtaMinCal;
 217     fEtaMaxCal = other.fEtaMaxCal;
 218     fPhiMinCal = other.fPhiMinCal;
 219     fPhiMaxCal = other.fPhiMaxCal;
 220     fHadCorr = other.fHadCorr;
 221     fHCorrection = other.fHCorrection;
 222     fNTracks = other.fNTracks;
 223     fNTracksCut = other.fNTracksCut;
 224     fOpt = other.fOpt;
 225     fDZ = other.fDZ;
 226     fDebug = other.fDebug;
 227     fReaderHeader = other.fReaderHeader;
 228     fMomentumArray = other.fMomentumArray;
 229     fUnitArray = other.fUnitArray;
 230     fRefArray = other.fRefArray;
 231     fProcId = other.fProcId;
 232     fTPCGrid = other.fTPCGrid;
 233     fEMCalGrid = other.fEMCalGrid;
 234     fGeom = other.fGeom;
 235     fESD = other.fESD;
 236     fGrid0 = other.fGrid0;
 237     fGrid1 = other.fGrid1;
 238     fGrid2 = other.fGrid2;
 239     fGrid3 = other.fGrid3;
 240     fGrid4 = other.fGrid4;
 241     fNphi = other.fNphi;
 242     fNeta = other.fNeta;
 243     fPhi2 = other.fPhi2;
 244     fEta2 = other.fEta2;
 245     fPhi = other.fPhi;
 246     fEta = other.fEta;
 247     fIndex = other.fIndex;
 248     fParams = other.fParams;
 249     fGrid = other.fGrid;
 250     fPhiMin = other.fPhiMin;
 251     fPhiMax = other.fPhiMax;
 252     fEtaMin = other.fEtaMin;
 253     fEtaMax = other.fEtaMax;
 254     fEtaBinInTPCAcc = other.fEtaBinInTPCAcc;
 255     fPhiBinInTPCAcc = other.fPhiBinInTPCAcc;
 256     fEtaBinInEMCalAcc = other.fEtaBinInEMCalAcc;
 257     fPhiBinInEMCalAcc = other.fPhiBinInEMCalAcc;
 258     fNbinPhi = other.fNbinPhi;
 259
 260     return (*this);
 261 }
 262
 263 //____________________________________________________________________________
 264 void AliJetFillUnitArrayTracks::InitParameters()
 265 {
 266   fHadCorr        = 0;             // For hadron correction
 267   fNumUnits = fGeom->GetNCells();  // Number of towers in EMCAL
 268   //fDebug = fReaderHeader->GetDebug();
 269
 270   fEtaMinCal = fGeom->GetArm1EtaMin();
 271   fEtaMaxCal = fGeom->GetArm1EtaMax();
 272   fPhiMinCal = fGeom->GetArm1PhiMin();
 273   fPhiMaxCal = fGeom->GetArm1PhiMax()-10.;
 274
 275   fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
 276                         fPhiMax,fEtaMin,fEtaMax);
 277   fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
 278                         fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
 279
 280   fEta   = fTPCGrid->GetArrayEta();
 281   fPhi   = fTPCGrid->GetArrayPhi();
 282   fIndex = fTPCGrid->GetIndexObject();
 283
 284   if(fDebug>20){
 285     for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
 286     for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
 287
 288     for(Int_t i=0; i<fNphi+1; i++)
 289       for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
 290           (*fIndex)(i,j) << endl; }
 291   }
 292   if(fDebug>1) printf("\n Parameters initiated ! \n");
 293 }
 294
 295 //_____________________________________________________________________________
 296 AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
 297 {
 298   // destructor
 299 }
 300
 301 //_____________________________________________________________________________
 302 void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
 303 {
 304   //
 305   // Main method.
 306   //
 307
 308   fDebug = fReaderHeader->GetDebug();
 309
 310   // Set parameters
 311   InitParameters();
 312
 313   // get number of tracks in event (for the loop)
 314   Int_t goodTrack = 0;
 315   Int_t nt = 0;
 316   Int_t nmax = 0;
 317   Float_t pt,eta,phi;
 318   //  Int_t sflag = 0;
 319   TVector3 p3;
 320
 321   if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl;
 322
 323   nt = fESD->GetNumberOfTracks();
 324   if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
 325
 326   // temporary storage of signal and pt cut flag
 327   Int_t* sflag  = new Int_t[nt];
 328   Int_t* cflag  = new Int_t[nt];
 329
 330   // get cuts set by user
 331   Float_t ptMin  = fReaderHeader->GetPtCut();
 332   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 333   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 334   Float_t phiMin = fReaderHeader->GetFiducialPhiMin();
 335   Float_t phiMax = fReaderHeader->GetFiducialPhiMax();
 336   fOpt = fReaderHeader->GetDetector();
 337   fDZ  = fReaderHeader->GetDZ();
 338
 339   Int_t nTracksEmcal      = 0;
 340   Int_t nTracksEmcalDZ    = 0;
 341   Int_t nTracksTpc        = 0;
 342   Int_t nTracksTpcOnly    = 0;
 343   Int_t nTracksEmcalCut   = 0;
 344   Int_t nTracksEmcalDZCut = 0;
 345   Int_t nTracksTpcCut     = 0;
 346   Int_t nTracksTpcOnlyCut = 0;
 347
 348   fGrid = fTPCGrid->GetGridType();
 349
 350   // Loop over tracks
 351   nmax = nt;
 352   for (Int_t it = 0; it < nmax; it++) {
 353     AliESDtrack *track;
 354     track = fESD->GetTrack(it);
 355     UInt_t status = track->GetStatus();
 356
 357     Double_t mom[3];
 358     track->GetPxPyPz(mom);
 359
 360     p3.SetXYZ(mom[0],mom[1],mom[2]);
 361     vector<Float_t> vtmp(3);
 362     vtmp[0] = mom[0]; vtmp[1] = mom[1]; vtmp[2] = mom[2];
 363     pt = p3.Pt();
 364
 365     Float_t mass = 0.;
 366     mass = track->GetMass();
 367
 368     if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
 369     if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
 370     if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
 371         && TMath::Abs(track->GetLabel()) > 10000)  continue;   // quality check
 372     if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
 373         && TMath::Abs(track->GetLabel()) < 10000)  continue;   // quality check
 374     eta = p3.Eta();
 375     phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
 376
 377     if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
 378     if ( (phi > phiMax) || (phi < phiMin)) continue;           // checking phi cut
 379
 380     sflag[goodTrack]=0;
 381     if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 382     cflag[goodTrack]=0;
 383     if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
 384
 385     if(fDebug>20) cout << "In Pythia: Track:" << it << ", eta: " << eta << ", phi: " << phi << ", pt: " << pt << endl;
 386
 387     if(fGrid==0)
 388       {
 389         // Only TPC filled from its grid in its total acceptance
 390
 391         Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
 392         Bool_t ok = kFALSE;
 393         Bool_t ref = kFALSE;
 394
 395         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
 396
 397         uArray->SetUnitTrackID(it);
 398
 399         Float_t unitEnergy = 0.;
 400         unitEnergy = uArray->GetUnitEnergy();
 401         // nTracksTpcOnly is necessary to count the number of candidate cells
 402         // but it doesn't give the real multiplicity -> it will be extracted
 403         // in the jet finder and for that, it is necessary to fill the Px,Py,Pz
 404         // information for each tracks stored in a given unitcell
 405         if(unitEnergy==0.){
 406           nTracksTpcOnly++;
 407           ok = kTRUE;
 408           ref = kTRUE;
 409         }
 410
 411         // Fill energy in TPC acceptance
 412         uArray->SetUnitEnergy(unitEnergy + pt);
 413         uArray->SetUnitPxPyPz(kFALSE,vtmp);
 414         uArray->SetUnitMass(mass);
 415
 416         // Pt cut flag
 417         if(uArray->GetUnitEnergy()<ptMin){
 418           uArray->SetUnitCutFlag(kPtSmaller);
 419         }
 420         else {
 421           uArray->SetUnitCutFlag(kPtHigher);
 422           if(ok) nTracksTpcOnlyCut++;
 423         }
 424
 425         // Signal flag
 426         if(sflag[goodTrack] == 1) {
 427           uArray->SetUnitSignalFlag(kGood);
 428           uArray->SetUnitSignalFlagC(kFALSE,kGood);
 429         } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
 430
 431         if(uArray->GetUnitEnergy()>0 && ref){
 432           if(!fProcId) {
 433             fProcId = kTRUE;
 434             // delete fRefArray;
 435             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 436           }
 437           fRefArray->Add(uArray);
 438         }
 439
 440       }
 441
 442     if(fGrid==1)
 443       {
 444         Int_t nElements = fTPCGrid->GetNEntries();
 445         // Fill track information in EMCAL acceptance
 446         if((eta >= fEtaMin && eta <= fEtaMax) &&
 447            (phi >= fPhiMin && phi <= fPhiMax))// &&
 448           {
 449             // Include dead-zones
 450             if(fDZ)
 451               {
 452                 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 453                 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 454                 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 455                 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 456                 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 457                 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 458                 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 459                 Int_t n0 = fGrid0->GetNEntries();
 460                 Int_t n1 = fGrid1->GetNEntries();
 461                 Int_t n2 = fGrid2->GetNEntries();
 462                 Int_t n3 = fGrid3->GetNEntries();
 463
 464                 if(phi >= phimin0 && phi <= phimax0){
 465                   Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
 466                   AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
 467
 468                   uArray0->SetUnitTrackID(it);
 469
 470                   Float_t uEnergy0 = uArray0->GetUnitEnergy();
 471                   Bool_t ok0 = kFALSE;
 472                   Bool_t ref0 = kFALSE;
 473                   if(uEnergy0==0.){
 474                     nTracksEmcalDZ++;
 475                     ok0 = kTRUE;
 476                     ref0 = kTRUE;
 477                   }
 478                   uArray0->SetUnitEnergy(uEnergy0+pt);
 479                   uArray0->SetUnitPxPyPz(kFALSE,vtmp);
 480                   if(uArray0->GetUnitEnergy()<ptMin)
 481                     uArray0->SetUnitCutFlag(kPtSmaller);
 482                   else {
 483                     uArray0->SetUnitCutFlag(kPtHigher);
 484                     if(ok0) nTracksEmcalDZCut++;
 485                   }
 486                   if(sflag[goodTrack] == 1) {
 487                     uArray0->SetUnitSignalFlag(kGood);
 488                     uArray0->SetUnitSignalFlagC(kFALSE,kGood);
 489                   } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
 490
 491                   if(uArray0->GetUnitEnergy()>0 && ref0){
 492                     if(!fProcId){
 493                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
 494                       fProcId = kTRUE;
 495                     }
 496                     fRefArray->Add(uArray0);
 497                   }
 498                 }
 499                 if(phi >= phimin1 && phi <= phimax1){
 500                   Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
 501                   AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
 502                   uArray1->SetUnitTrackID(it);
 503
 504                   Float_t uEnergy1 = uArray1->GetUnitEnergy();
 505                   Bool_t ok1 = kFALSE;
 506                   Bool_t ref1 = kFALSE;
 507                   if(uEnergy1==0.){
 508                     nTracksEmcalDZ++;
 509                     ok1 = kTRUE;
 510                     ref1 = kTRUE;
 511                   }
 512                   uArray1->SetUnitEnergy(uEnergy1+pt);
 513                   uArray1->SetUnitPxPyPz(kFALSE,vtmp);
 514                   if(uArray1->GetUnitEnergy()<ptMin)
 515                     uArray1->SetUnitCutFlag(kPtSmaller);
 516                   else {
 517                     uArray1->SetUnitCutFlag(kPtHigher);
 518                     if(ok1) nTracksEmcalDZCut++;
 519                   }
 520                   if(sflag[goodTrack] == 1) {
 521                     uArray1->SetUnitSignalFlag(kGood);
 522                     uArray1->SetUnitSignalFlagC(kFALSE,kGood);
 523                   } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
 524
 525                   if(uArray1->GetUnitEnergy()>0 && ref1){
 526                     if(!fProcId){
 527                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
 528                       fProcId = kTRUE;
 529                     }
 530                     fRefArray->Add(uArray1);
 531                   }
 532                 }
 533                 if(phi >= phimin2 && phi <= phimax2){
 534                   Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
 535                   AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
 536                   uArray2->SetUnitTrackID(it);
 537
 538                   Float_t uEnergy2 = uArray2->GetUnitEnergy();
 539                   Bool_t ok2 = kFALSE;
 540                   Bool_t ref2 = kFALSE;
 541                   if(uEnergy2==0.){
 542                     nTracksEmcalDZ++;
 543                     ok2 = kTRUE;
 544                     ref2 = kTRUE;
 545                   }
 546                   uArray2->SetUnitEnergy(uEnergy2+pt);
 547                   uArray2->SetUnitPxPyPz(kFALSE,vtmp);
 548                   if(uArray2->GetUnitEnergy()<ptMin)
 549                     uArray2->SetUnitCutFlag(kPtSmaller);
 550                   else {
 551                     uArray2->SetUnitCutFlag(kPtHigher);
 552                     if(ok2) nTracksEmcalDZCut++;
 553                   }
 554                   if(sflag[goodTrack] == 1) {
 555                     uArray2->SetUnitSignalFlag(kGood);
 556                     uArray2->SetUnitSignalFlagC(kFALSE,kGood);
 557                   } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
 558
 559                   if(uArray2->GetUnitEnergy()>0 && ref2){
 560                     if(!fProcId){
 561                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
 562                       fProcId = kTRUE;
 563                     }
 564                     fRefArray->Add(uArray2);
 565                   }
 566                 }
 567                 if(phi >= phimin3 && phi <= phimax3){
 568                   Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
 569                   AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
 570                   uArray3->SetUnitTrackID(it);
 571
 572                   Float_t uEnergy3 = uArray3->GetUnitEnergy();
 573                   Bool_t ok3 = kFALSE;
 574                   Bool_t ref3 = kFALSE;
 575                   if(uEnergy3==0.){
 576                     nTracksEmcalDZ++;
 577                     ok3 = kTRUE;
 578                     ref3 = kTRUE;
 579                   }
 580                   uArray3->SetUnitEnergy(uEnergy3+pt);
 581                   uArray3->SetUnitPxPyPz(kFALSE,vtmp);
 582                   if(uArray3->GetUnitEnergy()<ptMin)
 583                     uArray3->SetUnitCutFlag(kPtSmaller);
 584                   else {
 585                     uArray3->SetUnitCutFlag(kPtHigher);
 586                     if(ok3) nTracksEmcalDZCut++;
 587                   }
 588                   if(sflag[goodTrack] == 1) {
 589                     uArray3->SetUnitSignalFlag(kGood);
 590                     uArray3->SetUnitSignalFlagC(kFALSE,kGood);
 591                   } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
 592
 593                   if(uArray3->GetUnitEnergy()>0 && ref3){
 594                     if(!fProcId){
 595                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
 596                       fProcId = kTRUE;
 597                     }
 598                     fRefArray->Add(uArray3);
 599                   }
 600                 }
 601                 if(phi >= phimin4 && phi <= phimax4){
 602                   Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
 603                   AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
 604                   uArray4->SetUnitTrackID(it);
 605
 606                   Float_t uEnergy4 = uArray4->GetUnitEnergy();
 607                   Bool_t ok4 = kFALSE;
 608                   Bool_t ref4 = kFALSE;
 609                   if(uEnergy4==0.){
 610                     nTracksEmcalDZ++;
 611                     ok4 = kTRUE;
 612                     ref4 = kTRUE;
 613                   }
 614                   uArray4->SetUnitEnergy(uEnergy4+pt);
 615                   uArray4->SetUnitPxPyPz(kFALSE,vtmp);
 616                   if(uArray4->GetUnitEnergy()<ptMin)
 617                     uArray4->SetUnitCutFlag(kPtSmaller);
 618                   else {
 619                     uArray4->SetUnitCutFlag(kPtHigher);
 620                     if(ok4) nTracksEmcalDZCut++;
 621                   }
 622                   if(sflag[goodTrack] == 1) {
 623                     uArray4->SetUnitSignalFlag(kGood);
 624                     uArray4->SetUnitSignalFlagC(kFALSE,kGood);
 625                   } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
 626
 627                   if(uArray4->GetUnitEnergy()>0 && ref4){
 628                     if(!fProcId){
 629                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
 630                       fProcId = kTRUE;
 631                     }
 632                     fRefArray->Add(uArray4);
 633                   }
 634                 }
 635               } // end fDZ
 636
 637             Int_t towerID = 0;
 638             fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
 639             if(towerID==-1) continue;
 640
 641             AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
 642             uArray->SetUnitTrackID(it);
 643
 644             Float_t 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             // Do Hadron Correction
 654             // Parametrization to be added
 655             if (fHCorrection != 0)
 656               {
 657                 Float_t   hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta);
 658                 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
 659
 660               } //end Hadron Correction loop
 661
 662             uArray->SetUnitEnergy(unitEnergy + pt);
 663             uArray->SetUnitPxPyPz(kFALSE,vtmp);
 664             // Put a pt cut flag
 665             if(uArray->GetUnitEnergy()<ptMin){
 666               uArray->SetUnitCutFlag(kPtSmaller);
 667             }
 668             else {
 669               uArray->SetUnitCutFlag(kPtHigher);
 670               if(ok) nTracksEmcalCut++;
 671             }
 672
 673             // Signal flag
 674             if(sflag[goodTrack] == 1) {
 675               uArray->SetUnitSignalFlag(kGood);
 676               uArray->SetUnitSignalFlagC(kFALSE,kGood);
 677             } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
 678
 679             if(uArray->GetUnitEnergy()>0 && ref){
 680               if(!fProcId){
 681                 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 682                 fProcId = kTRUE;
 683               }
 684               fRefArray->Add(uArray);
 685             }
 686           } // end loop on EMCal acceptance cut + tracks quality
 687         else{
 688           // Outside EMCal acceptance
 689
 690           Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
 691
 692           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
 693           uArray->SetUnitTrackID(it);
 694
 695           Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
 696           Bool_t okout = kFALSE;
 697           Bool_t refout = kFALSE;
 698           if(unitEnergy2==0.){
 699             nTracksTpc++;
 700             okout=kTRUE;
 701             refout=kTRUE;
 702           }
 703           // Fill energy outside emcal acceptance
 704           uArray->SetUnitEnergy(unitEnergy2 + pt);
 705           uArray->SetUnitPxPyPz(kFALSE,vtmp);
 706
 707           // Pt cut flag
 708           if(uArray->GetUnitEnergy()<ptMin){
 709             uArray->SetUnitCutFlag(kPtSmaller);
 710           }
 711           else {
 712             uArray->SetUnitCutFlag(kPtHigher);
 713             if(okout) nTracksTpcCut++;
 714           }
 715           // Signal flag
 716           if(sflag[goodTrack] == 1) {
 717             uArray->SetUnitSignalFlag(kGood);
 718             uArray->SetUnitSignalFlagC(kFALSE,kGood);
 719           } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
 720
 721           if(uArray->GetUnitEnergy()>0 && refout){
 722             if(!fProcId){
 723               new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 724               fProcId = kTRUE;
 725             }
 726             fRefArray->Add(uArray);
 727           }
 728         }
 729       } // end fGrid==1
 730
 731     goodTrack++;
 732
 733   } // end loop on entries (tpc tracks)
 734
 735   if(fDebug>0)
 736     {
 737       cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
 738       cout << "goodTracks: " << goodTrack << endl;
 739     }
 740
 741   delete[] sflag;
 742   delete[] cflag;
 743
 744   if(fGrid==0) {
 745     fNTracks = nTracksTpcOnly;
 746     fNTracksCut = nTracksTpcOnlyCut;
 747     if(fDebug>10){
 748       cout << "fNTracks : " << fNTracks << endl;
 749       cout << "fNTracksCut : " << fNTracksCut << endl;
 750     }
 751   }
 752   if(fGrid==1) {
 753     fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
 754     fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
 755     if(fDebug>10){
 756       cout << "fNTracks : " << fNTracks << endl;
 757       cout << "fNTracksCut : " << fNTracksCut << endl;
 758     }
 759   }
 760
 761 }
 762
 763 //__________________________________________________________
 764 void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
 765 {
 766   for(Int_t j=0; j<fNphi+1; j++)
 767     {
 768       for(Int_t i=0; i<fNeta+1; i++)
 769         {
 770
 771           // TPC grid only
 772           //-------------------------------------
 773           if(fGrid==0) {
 774             if(j*(fNeta+1)+i == index) {
 775               eta = fEta2->At(i);
 776               phi = fPhi2->At(j);
 777             }
 778           }
 779
 780           // TPC-EMCAL grid
 781           //-------------------------------------
 782           Int_t ii = 0;
 783           if(i==0) ii = 0;
 784           if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
 785           if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
 786
 787           if(fGrid==1) {
 788             if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
 789               eta = fEta2->At(i);
 790               phi = fPhi2->At(j);
 791             }
 792
 793             if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
 794                ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
 795               if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
 796               else eta = fEta2->At(i);
 797               phi = fPhi2->At(j);
 798             }
 799
 800             if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
 801               eta = fEta2->At(i);
 802               phi = fPhi2->At(j);
 803             }
 804           }
 805         }
 806     }
 807 }
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819