]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - JETAN/AliJetFillUnitArrayTracks.cxx
Updates needed for full jet reconstruction (charged + emcal) [Magali Estienne]
[u/mrichter/AliRoot.git] / JETAN / AliJetFillUnitArrayTracks.cxx
index 18ab76657fd70d4844e8f99eb92720523b5e3111..b96ad08d45d0d611af52f1c1f57882fd4a3cd6fa 100644 (file)
  **************************************************************************/
 
 
-//----------------------------------------------------------------------
+//======================================================================
+// ***July 2006
 // Fill Unit Array class 
-// Class used by AliJetESDReader to fill a UnitArray from the information
-// extracted from the particle tracks
+// Class used by AliJetESDReader to fill a UnitArray from the information extracted 
+// from the particle tracks
 // Author: magali.estienne@ires.in2p3.fr
-//----------------------------------------------------------------------
+//======================================================================
 
 
 // --- Standard library ---
@@ -36,6 +37,7 @@
 #include <TArrayD.h>
 #include <TMath.h>
 #include <TClonesArray.h>
+#include <TProcessID.h>
 
 // --- AliRoot header files ---
 #include "AliJetFinder.h"
@@ -43,7 +45,6 @@
 #include "AliJetReader.h"
 #include "AliJetESDReader.h"
 #include "AliJetESDReaderHeader.h"
-//#include "AliESD.h"
 #include "AliESDEvent.h"
 #include "AliJetDummyGeo.h"
 #include "AliJetUnitArray.h"
@@ -72,6 +73,7 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
     fMomentumArray(0x0),
     fUnitArray(0x0),
     fRefArray(0x0),
+    fProcId(kFALSE),
     fTPCGrid(0x0),
     fEMCalGrid(0x0),
     fGeom(0x0),
@@ -104,7 +106,7 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
 }
 
 //_____________________________________________________________________________
-AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/)
+AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* esd)
   : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
     fNumUnits(0),
     fEtaMinCal(0),
@@ -122,15 +124,16 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/)
     fMomentumArray(0x0),
     fUnitArray(0x0),
     fRefArray(0x0),
+    fProcId(kFALSE),
     fTPCGrid(0x0),
     fEMCalGrid(0x0),
     fGeom(0x0),
-    fESD(0x0),
+    fESD(esd),
     fGrid0(0x0),
     fGrid1(0x0),
     fGrid2(0x0),
     fGrid3(0x0),
-    fGrid4(0x0),
+    fGrid4(0x0),// first: loop over tracks in ESD 
     fNphi(0),
     fNeta(0),
     fPhi2(0x0),
@@ -150,21 +153,124 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/)
     fPhiBinInEMCalAcc(0),
     fNbinPhi(0)
 {
-  // constructor
+  // constructor 2
+}
+
+//_____________________________________________________________________________
+AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(const AliJetFillUnitArrayTracks &det)
+  : TTask(det),//"AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
+    fNumUnits(det.fNumUnits),
+    fEtaMinCal(det.fEtaMinCal),
+    fEtaMaxCal(det.fEtaMaxCal),
+    fPhiMinCal(det.fPhiMinCal),
+    fPhiMaxCal(det.fPhiMaxCal),
+    fHadCorr(det.fHadCorr),
+    fHCorrection(det.fHCorrection),
+    fNTracks(det.fNTracks),
+    fNTracksCut(det.fNTracksCut),
+    fOpt(det.fOpt),
+    fDZ(det.fDZ),
+    fDebug(det.fDebug),
+    fReaderHeader(det.fReaderHeader),
+    fMomentumArray(det.fMomentumArray),
+    fUnitArray(det.fUnitArray),
+    fRefArray(det.fRefArray),
+    fProcId(det.fProcId),
+    fTPCGrid(det.fTPCGrid),
+    fEMCalGrid(det.fEMCalGrid),
+    fGeom(det.fGeom),
+    fESD(det.fESD),
+    fGrid0(det.fGrid0),
+    fGrid1(det.fGrid1),
+    fGrid2(det.fGrid2),
+    fGrid3(det.fGrid3),
+    fGrid4(det.fGrid4),
+    fNphi(det.fNphi),
+    fNeta(det.fNeta),
+    fPhi2(det.fPhi2),
+    fEta2(det.fEta2),
+    fPhi(det.fPhi),
+    fEta(det.fEta),
+    fIndex(det.fIndex),
+    fParams(det.fParams),
+    fGrid(det.fGrid),
+    fPhiMin(det.fPhiMin),
+    fPhiMax(det.fPhiMax),
+    fEtaMin(det.fEtaMin),
+    fEtaMax(det.fEtaMax),
+    fEtaBinInTPCAcc(det.fEtaBinInTPCAcc),
+    fPhiBinInTPCAcc(det.fPhiBinInTPCAcc),
+    fEtaBinInEMCalAcc(det.fEtaBinInEMCalAcc),
+    fPhiBinInEMCalAcc(det.fPhiBinInEMCalAcc),
+    fNbinPhi(det.fNbinPhi)
+{
+  // Copy constructor 
+}
+
+//------------------------------------------------------------------
+AliJetFillUnitArrayTracks& AliJetFillUnitArrayTracks::operator=(const AliJetFillUnitArrayTracks& other)
+{
+  // Assignment
+
+    fNumUnits = other.fNumUnits;
+    fEtaMinCal = other.fEtaMinCal;
+    fEtaMaxCal = other.fEtaMaxCal;
+    fPhiMinCal = other.fPhiMinCal;
+    fPhiMaxCal = other.fPhiMaxCal;
+    fHadCorr = other.fHadCorr;
+    fHCorrection = other.fHCorrection;
+    fNTracks = other.fNTracks;
+    fNTracksCut = other.fNTracksCut;
+    fOpt = other.fOpt;
+    fDZ = other.fDZ;
+    fDebug = other.fDebug;
+    fReaderHeader = other.fReaderHeader;
+    fMomentumArray = other.fMomentumArray;
+    fUnitArray = other.fUnitArray;
+    fRefArray = other.fRefArray;
+    fProcId = other.fProcId;
+    fTPCGrid = other.fTPCGrid;
+    fEMCalGrid = other.fEMCalGrid;
+    fGeom = other.fGeom;
+    fESD = other.fESD;
+    fGrid0 = other.fGrid0;
+    fGrid1 = other.fGrid1;
+    fGrid2 = other.fGrid2;
+    fGrid3 = other.fGrid3;
+    fGrid4 = other.fGrid4;
+    fNphi = other.fNphi;
+    fNeta = other.fNeta;
+    fPhi2 = other.fPhi2;
+    fEta2 = other.fEta2;
+    fPhi = other.fPhi;
+    fEta = other.fEta;
+    fIndex = other.fIndex;
+    fParams = other.fParams;
+    fGrid = other.fGrid;
+    fPhiMin = other.fPhiMin;
+    fPhiMax = other.fPhiMax;
+    fEtaMin = other.fEtaMin;
+    fEtaMax = other.fEtaMax;
+    fEtaBinInTPCAcc = other.fEtaBinInTPCAcc;
+    fPhiBinInTPCAcc = other.fPhiBinInTPCAcc;
+    fEtaBinInEMCalAcc = other.fEtaBinInEMCalAcc;
+    fPhiBinInEMCalAcc = other.fPhiBinInEMCalAcc;
+    fNbinPhi = other.fNbinPhi;
+
+    return (*this);
 }
 
 //____________________________________________________________________________
 void AliJetFillUnitArrayTracks::InitParameters()
 {
-  //  fHCorrection    = 0;     // For hadron correction
-  fHadCorr        = 0;     // For hadron correction
-  fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
-  fDebug = fReaderHeader->GetDebug();
+  fHadCorr        = 0;             // For hadron correction
+  fNumUnits = fGeom->GetNCells();  // Number of towers in EMCAL
+  //fDebug = fReaderHeader->GetDebug();
 
   fEtaMinCal = fGeom->GetArm1EtaMin();
   fEtaMaxCal = fGeom->GetArm1EtaMax();
   fPhiMinCal = fGeom->GetArm1PhiMin();
-  fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
+  fPhiMaxCal = fGeom->GetArm1PhiMax()-10.;
 
   fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, 
                        fPhiMax,fEtaMin,fEtaMax);
@@ -200,19 +306,24 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
   //
 
   fDebug = fReaderHeader->GetDebug();
-
+  
   // Set parameters
   InitParameters();
 
   // get number of tracks in event (for the loop)
   Int_t goodTrack = 0;
   Int_t nt = 0;
-  Float_t pt, eta,phi;
+  Int_t nmax = 0;
+  Float_t pt,eta,phi;
+  //  Int_t sflag = 0;
   TVector3 p3;
+
+  if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl;
+
   nt = fESD->GetNumberOfTracks();
   if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
-   
- // temporary storage of signal and pt cut flag
 // temporary storage of signal and pt cut flag
   Int_t* sflag  = new Int_t[nt];
   Int_t* cflag  = new Int_t[nt];
 
@@ -220,6 +331,8 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
   Float_t ptMin  = fReaderHeader->GetPtCut();
   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
+  Float_t phiMin = fReaderHeader->GetFiducialPhiMin();
+  Float_t phiMax = fReaderHeader->GetFiducialPhiMax();  
   fOpt = fReaderHeader->GetDetector();
   fDZ  = fReaderHeader->GetDZ();
 
@@ -234,52 +347,71 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
 
   fGrid = fTPCGrid->GetGridType();
 
-
-  //loop over tracks
-  for (Int_t it = 0; it < nt; it++) {
-    AliESDtrack *track = fESD->GetTrack(it);
+  // Loop over tracks
+  nmax = nt;  
+  for (Int_t it = 0; it < nmax; it++) {
+    AliESDtrack *track;
+    track = fESD->GetTrack(it);
     UInt_t status = track->GetStatus();
     
     Double_t mom[3];
     track->GetPxPyPz(mom);
+    
     p3.SetXYZ(mom[0],mom[1],mom[2]);
+    vector<Float_t> vtmp(3);
+    vtmp[0] = mom[0]; vtmp[1] = mom[1]; vtmp[2] = mom[2];
     pt = p3.Pt();
+    
     Float_t mass = 0.;
     mass = track->GetMass();
-    
-    if (((status & AliESDtrack::kITSrefit) == 0) ||
-        ((status & AliESDtrack::kTPCrefit) == 0)) continue;    // quality check
+
+    if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
+    if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
     if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() 
        && TMath::Abs(track->GetLabel()) > 10000)  continue;   // quality check
     if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() 
        && TMath::Abs(track->GetLabel()) < 10000)  continue;   // quality check
     eta = p3.Eta();
     phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
-    
+       
     if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
+    if ( (phi > phiMax) || (phi < phiMin)) continue;           // checking phi cut
+      
+    sflag[goodTrack]=0;
+    if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+    cflag[goodTrack]=0;
+    if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
+
+    if(fDebug>20) cout << "In Pythia: Track:" << it << ", eta: " << eta << ", phi: " << phi << ", pt: " << pt << endl;
 
-    // sflag -> not yet implemented !!!!
-    
     if(fGrid==0)
       {
        // Only TPC filled from its grid in its total acceptance
        
        Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
        Bool_t ok = kFALSE;
-
+       Bool_t ref = kFALSE;
+       
        AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
+       
        uArray->SetUnitTrackID(it);
-
+       
        Float_t unitEnergy = 0.;
        unitEnergy = uArray->GetUnitEnergy();
+       // nTracksTpcOnly is necessary to count the number of candidate cells
+       // but it doesn't give the real multiplicity -> it will be extracted 
+       // in the jet finder and for that, it is necessary to fill the Px,Py,Pz
+       // information for each tracks stored in a given unitcell 
        if(unitEnergy==0.){
          nTracksTpcOnly++;
          ok = kTRUE;
+         ref = kTRUE;
        }
+
        // Fill energy in TPC acceptance
        uArray->SetUnitEnergy(unitEnergy + pt);
-       uArray->SetUnitPxPyPz(mom);
-       uArray->SetUnitMass(mass);
+       uArray->SetUnitPxPyPz(kFALSE,vtmp);
+       uArray->SetUnitMass(mass);
 
        // Pt cut flag
        if(uArray->GetUnitEnergy()<ptMin){
@@ -290,31 +422,30 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
          if(ok) nTracksTpcOnlyCut++;
        }
 
-       // Detector flag
-       if(unitEnergy>0) {
-         uArray->SetUnitDetectorFlag(kAll);
-       }
-       if(uArray->GetUnitEnergy()>0){
+       // Signal flag
+       if(sflag[goodTrack] == 1) { 
+         uArray->SetUnitSignalFlag(kGood);
+         uArray->SetUnitSignalFlagC(kFALSE,kGood);
+       } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
+
+       if(uArray->GetUnitEnergy()>0 && ref){
+         if(!fProcId) {
+           fProcId = kTRUE;
+           // delete fRefArray;
+           new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+         }
          fRefArray->Add(uArray);
        }
-
-       sflag[goodTrack]=0;
-       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
-       cflag[goodTrack]=0;
-       if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
-       goodTrack++;
-
+       
       }
     
     if(fGrid==1)
       {
        Int_t nElements = fTPCGrid->GetNEntries();
-
        // Fill track information in EMCAL acceptance
        if((eta >= fEtaMin && eta <= fEtaMax) &&
           (phi >= fPhiMin && phi <= fPhiMax))// &&
          {
-
            // Include dead-zones
            if(fDZ)
              {
@@ -329,173 +460,249 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
                Int_t n1 = fGrid1->GetNEntries();
                Int_t n2 = fGrid2->GetNEntries();
                Int_t n3 = fGrid3->GetNEntries();
-
+               
                if(phi >= phimin0 && phi <= phimax0){
                  Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
                  AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
+                 
                  uArray0->SetUnitTrackID(it);
+                 
                  Float_t uEnergy0 = uArray0->GetUnitEnergy();
                  Bool_t ok0 = kFALSE;
+                 Bool_t ref0 = kFALSE;
                  if(uEnergy0==0.){
                    nTracksEmcalDZ++;
                    ok0 = kTRUE;
+                   ref0 = kTRUE;
                  }
                  uArray0->SetUnitEnergy(uEnergy0+pt);
+                 uArray0->SetUnitPxPyPz(kFALSE,vtmp);
                  if(uArray0->GetUnitEnergy()<ptMin)
                    uArray0->SetUnitCutFlag(kPtSmaller);
                  else {
                    uArray0->SetUnitCutFlag(kPtHigher);
                    if(ok0) nTracksEmcalDZCut++;
                  }
-                 if(uArray0->GetUnitEnergy()>0)
+                 if(sflag[goodTrack] == 1) {
+                   uArray0->SetUnitSignalFlag(kGood);
+                   uArray0->SetUnitSignalFlagC(kFALSE,kGood);
+                 } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
+                 
+                 if(uArray0->GetUnitEnergy()>0 && ref0){
+                   if(!fProcId){
+                     new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
+                     fProcId = kTRUE;
+                   }
                    fRefArray->Add(uArray0);
+                 }
                }
                if(phi >= phimin1 && phi <= phimax1){
                  Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
                  AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
                  uArray1->SetUnitTrackID(it);
+
                  Float_t uEnergy1 = uArray1->GetUnitEnergy();
                  Bool_t ok1 = kFALSE;
+                 Bool_t ref1 = kFALSE;
                  if(uEnergy1==0.){
                    nTracksEmcalDZ++;
                    ok1 = kTRUE;
+                   ref1 = kTRUE;
                  }
                  uArray1->SetUnitEnergy(uEnergy1+pt);
+                 uArray1->SetUnitPxPyPz(kFALSE,vtmp);
                  if(uArray1->GetUnitEnergy()<ptMin)
                    uArray1->SetUnitCutFlag(kPtSmaller);
                  else {
                    uArray1->SetUnitCutFlag(kPtHigher);
                    if(ok1) nTracksEmcalDZCut++;
                  }
-                 if(uArray1->GetUnitEnergy()>0)
+                 if(sflag[goodTrack] == 1) {
+                   uArray1->SetUnitSignalFlag(kGood);
+                   uArray1->SetUnitSignalFlagC(kFALSE,kGood);
+                 } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
+                 
+                 if(uArray1->GetUnitEnergy()>0 && ref1){
+                   if(!fProcId){
+                     new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
+                     fProcId = kTRUE;
+                   }
                    fRefArray->Add(uArray1);
+                 }
                }
                if(phi >= phimin2 && phi <= phimax2){
                  Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
                  AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
                  uArray2->SetUnitTrackID(it);
+
                  Float_t uEnergy2 = uArray2->GetUnitEnergy();
                  Bool_t ok2 = kFALSE;
+                 Bool_t ref2 = kFALSE;
                  if(uEnergy2==0.){
                    nTracksEmcalDZ++;
                    ok2 = kTRUE;
+                   ref2 = kTRUE;
                  }
                  uArray2->SetUnitEnergy(uEnergy2+pt);
+                 uArray2->SetUnitPxPyPz(kFALSE,vtmp);
                  if(uArray2->GetUnitEnergy()<ptMin)
                    uArray2->SetUnitCutFlag(kPtSmaller);
                  else {
                    uArray2->SetUnitCutFlag(kPtHigher);
                    if(ok2) nTracksEmcalDZCut++;
                  }
-                 if(uArray2->GetUnitEnergy()>0)
+                 if(sflag[goodTrack] == 1) {
+                   uArray2->SetUnitSignalFlag(kGood);
+                   uArray2->SetUnitSignalFlagC(kFALSE,kGood);
+                 } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
+                 
+                 if(uArray2->GetUnitEnergy()>0 && ref2){
+                   if(!fProcId){
+                     new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
+                     fProcId = kTRUE;
+                   }
                    fRefArray->Add(uArray2);
+                 }
                }
                if(phi >= phimin3 && phi <= phimax3){
                  Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
                  AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
                  uArray3->SetUnitTrackID(it);
+
                  Float_t uEnergy3 = uArray3->GetUnitEnergy();
                  Bool_t ok3 = kFALSE;
+                 Bool_t ref3 = kFALSE;
                  if(uEnergy3==0.){
                    nTracksEmcalDZ++;
                    ok3 = kTRUE;
+                   ref3 = kTRUE;
                  }
                  uArray3->SetUnitEnergy(uEnergy3+pt);
+                 uArray3->SetUnitPxPyPz(kFALSE,vtmp);
                  if(uArray3->GetUnitEnergy()<ptMin)
                    uArray3->SetUnitCutFlag(kPtSmaller);
                  else {
                    uArray3->SetUnitCutFlag(kPtHigher);
                    if(ok3) nTracksEmcalDZCut++;
                  }
-                 if(uArray3->GetUnitEnergy()>0)
+                 if(sflag[goodTrack] == 1) {
+                   uArray3->SetUnitSignalFlag(kGood);
+                   uArray3->SetUnitSignalFlagC(kFALSE,kGood);
+                 } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
+                 
+                 if(uArray3->GetUnitEnergy()>0 && ref3){
+                   if(!fProcId){
+                     new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
+                     fProcId = kTRUE;
+                   }
                    fRefArray->Add(uArray3);
+                 }
                }
                if(phi >= phimin4 && phi <= phimax4){
                  Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
                  AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
                  uArray4->SetUnitTrackID(it);
+                 
                  Float_t uEnergy4 = uArray4->GetUnitEnergy();
                  Bool_t ok4 = kFALSE;
+                 Bool_t ref4 = kFALSE;
                  if(uEnergy4==0.){
                    nTracksEmcalDZ++;
                    ok4 = kTRUE;
+                   ref4 = kTRUE;
                  }
                  uArray4->SetUnitEnergy(uEnergy4+pt);
+                 uArray4->SetUnitPxPyPz(kFALSE,vtmp);
                  if(uArray4->GetUnitEnergy()<ptMin)
                    uArray4->SetUnitCutFlag(kPtSmaller);
                  else {
                    uArray4->SetUnitCutFlag(kPtHigher);
                    if(ok4) nTracksEmcalDZCut++;
                  }
-                 if(uArray4->GetUnitEnergy()>0)
+                 if(sflag[goodTrack] == 1) {
+                   uArray4->SetUnitSignalFlag(kGood);
+                   uArray4->SetUnitSignalFlagC(kFALSE,kGood);
+                 } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
+                 
+                 if(uArray4->GetUnitEnergy()>0 && ref4){
+                   if(!fProcId){
+                     new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
+                     fProcId = kTRUE;
+                   }
                    fRefArray->Add(uArray4);
+                 }
                }
              } // end fDZ
            
            Int_t towerID = 0;
            fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
-
            if(towerID==-1) continue;
-
+             
            AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
            uArray->SetUnitTrackID(it);
+
            Float_t unitEnergy = uArray->GetUnitEnergy(); 
            Bool_t ok = kFALSE;
+           Bool_t ref = kFALSE;
            if(unitEnergy==0.){
              nTracksEmcal++;
              ok=kTRUE;
+             ref=kTRUE;
            }
 
            // Do Hadron Correction
            // Parametrization to be added
            if (fHCorrection != 0) 
              { 
-               //            Float_t   hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
                Float_t   hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta);
                unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
                
              } //end Hadron Correction loop
-
+           
            uArray->SetUnitEnergy(unitEnergy + pt);
-
+           uArray->SetUnitPxPyPz(kFALSE,vtmp);
            // Put a pt cut flag
-           if(uArray->GetUnitEnergy()<ptMin)
+           if(uArray->GetUnitEnergy()<ptMin){
              uArray->SetUnitCutFlag(kPtSmaller);
+           }
            else {
              uArray->SetUnitCutFlag(kPtHigher);
              if(ok) nTracksEmcalCut++;
            }
-           // Detector flag
-           if(unitEnergy > 0)
-             uArray->SetUnitDetectorFlag(kAll);
            
-           if(uArray->GetUnitEnergy()>0)
-             fRefArray->Add(uArray);
-               
-           sflag[goodTrack]=0;
-           if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
-           cflag[goodTrack]=0;
-           if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
-           goodTrack++;
+           // Signal flag
+           if(sflag[goodTrack] == 1) { 
+             uArray->SetUnitSignalFlag(kGood);
+             uArray->SetUnitSignalFlagC(kFALSE,kGood);
+           } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
            
+           if(uArray->GetUnitEnergy()>0 && ref){
+             if(!fProcId){
+               new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+               fProcId = kTRUE;
+             }
+             fRefArray->Add(uArray);
+           }
          } // end loop on EMCal acceptance cut + tracks quality
        else{ 
          // Outside EMCal acceptance
          
-         //    Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
          Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
-
+         
          AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
          uArray->SetUnitTrackID(it);
 
          Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
-         Bool_t ok2 = kFALSE;
+         Bool_t okout = kFALSE;
+         Bool_t refout = kFALSE;
          if(unitEnergy2==0.){
            nTracksTpc++;
-           ok2=kTRUE;
+           okout=kTRUE;
+           refout=kTRUE;
          }
          // Fill energy outside emcal acceptance
          uArray->SetUnitEnergy(unitEnergy2 + pt);
+         uArray->SetUnitPxPyPz(kFALSE,vtmp);
          
          // Pt cut flag
          if(uArray->GetUnitEnergy()<ptMin){
@@ -503,33 +710,37 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
          }
          else {
            uArray->SetUnitCutFlag(kPtHigher);
-           if(ok2) nTracksTpcCut++;
+           if(okout) nTracksTpcCut++;
          }
-         // Detector flag
-         if(unitEnergy2 > 0)
-           uArray->SetUnitDetectorFlag(kTpc);
-         if(uArray->GetUnitEnergy()>0)
-           fRefArray->Add(uArray);
-
-         sflag[goodTrack]=0;
-         if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
-         cflag[goodTrack]=0;
-         if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
-         goodTrack++;
+         // Signal flag
+         if(sflag[goodTrack] == 1) {
+           uArray->SetUnitSignalFlag(kGood);
+           uArray->SetUnitSignalFlagC(kFALSE,kGood);
+         } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
          
+         if(uArray->GetUnitEnergy()>0 && refout){
+           if(!fProcId){
+             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+             fProcId = kTRUE;
+           }
+           fRefArray->Add(uArray);
+         }
        }
       } // end fGrid==1
+
+    goodTrack++;
+
   } // end loop on entries (tpc tracks)
 
-//   // set the signal flags
-//   fSignalFlag.Set(goodTrack,sflag);
-//   fCutFlag.Set(goodTrack,cflag);
+  if(fDebug>0) 
+    {
+      cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
+      cout << "goodTracks: " << goodTrack << endl;
+    } 
 
-//   delete sflag;
-//   delete cflag;
+  delete[] sflag;
+  delete[] cflag;
 
-  //    } // end loop on entries (tpc tracks)
-      
   if(fGrid==0) {
     fNTracks = nTracksTpcOnly;
     fNTracksCut = nTracksTpcOnlyCut;
@@ -546,54 +757,53 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
       cout << "fNTracksCut : " << fNTracksCut << endl;
     }
   }  
-
+  
 }
 
 //__________________________________________________________
 void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
 {
-//
-// Obtain the index from eta and phi
-//
-  for(Int_t j=0; j<fNphi+1; j++) {
-    for(Int_t i=0; i<fNeta+1; i++) {
-
-      // TPC grid only 
-      //-------------------------------------
-      if(fGrid==0) {   
-       if(j*(fNeta+1)+i == index) {
-         eta = fEta2->At(i); 
-         phi = fPhi2->At(j);
-       }
-      }
-
-      // TPC-EMCAL grid
-      //-------------------------------------
-      Int_t ii = 0;
-      if(i==0) ii = 0;
-      if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i; 
-      if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
-
-      if(fGrid==1) {
-       if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
-         eta = fEta2->At(i);
-         phi = fPhi2->At(j);
-       }  
-
-       if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) && 
-          ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
-         if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
-         else eta = fEta2->At(i);
-         phi = fPhi2->At(j);
-       }
-
-       if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
-         eta = fEta2->At(i);
-         phi = fPhi2->At(j);
+  for(Int_t j=0; j<fNphi+1; j++) 
+    {
+      for(Int_t i=0; i<fNeta+1; i++) 
+       {
+
+         // TPC grid only 
+         //-------------------------------------
+         if(fGrid==0) {        
+           if(j*(fNeta+1)+i == index) {
+             eta = fEta2->At(i); 
+             phi = fPhi2->At(j);
+           }
+         }
+         
+         // TPC-EMCAL grid
+         //-------------------------------------
+         Int_t ii = 0;
+         if(i==0) ii = 0;
+         if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i; 
+         if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
+         
+         if(fGrid==1) {
+           if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
+             eta = fEta2->At(i);
+             phi = fPhi2->At(j);
+           }  
+           
+           if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) && 
+              ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
+             if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
+             else eta = fEta2->At(i);
+             phi = fPhi2->At(j);
+           }
+           
+           if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
+             eta = fEta2->At(i);
+             phi = fPhi2->At(j);
+           }
+         }
        }
-      }
     }
-  }
 }