[u/mrichter/AliRoot.git] / JETAN / AliJetESDReader.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//------------------------------------------------------------------------- 
// Jet ESD Reader 
// ESD reader for jet analysis
// Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
//          Magali Estienne <magali.estienne@subatech.in2p3.fr>
//
// **February 2011
// implemented standard geometry (AliEMCALGeometry) instead of dummy one (AliJetDummyGeo)
// moved geometry definition in AliJetReader
// marco.bregant@subatech.in2p3.fr
//------------------------------------------------------------------------- 

// --- Standard library ---
#include <Riostream.h>

// --- ROOT system ---
#include <TSystem.h>
#include <TString.h>
#include <TStopwatch.h>
#include <TLorentzVector.h>
#include <TVector3.h>
#include "TTask.h"
#include "TTree.h"
#include "TFile.h"
#include <TGeoManager.h>
#include <assert.h>
#include <TRefArray.h>
#include <TMath.h>
#include <TProcessID.h>
#include <TRandom3.h>

// --- AliRoot header files ---
#include "AliJetESDReader.h"
#include "AliJetESDReaderHeader.h"
#include "AliESDEvent.h"
#include "AliVEvent.h"
#include "AliESD.h"
#include "AliESDtrack.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALEMCGeometry.h"
#include "AliJetESDFillUnitArrayTracks.h"
#include "AliJetESDFillUnitArrayEMCalDigits.h"
#include "AliJetHadronCorrectionv1.h"
#include "AliJetUnitArray.h"
#include "AliAnalysisTask.h"
#include "AliOADBContainer.h"

using std::cout;
using std::endl;
ClassImp(AliJetESDReader)

AliJetESDReader::AliJetESDReader():
  AliJetReader(), 
  fHadCorr(0x0),
  fTpcGrid(0x0),
  fEmcalGrid(0x0),
  fGrid0(0),
  fGrid1(0),
  fGrid2(0),
  fGrid3(0),
  fGrid4(0),
  fApplyElectronCorrection(kFALSE),
  fApplyMIPCorrection(kTRUE),
  fApplyFractionHadronicCorrection(kFALSE),
  fFractionHadronicCorrection(0.3),
  fNumUnits(0),
  fDebug(0),
  fMass(0),
  fSign(0),
  fNIn(0),
  fOpt(0),
  fDZ(0),
  fNeta(0),
  fNphi(0),
  fRefArray(0x0),
  fProcId(kFALSE)
{
  // Constructor 
}

//____________________________________________________________________________
AliJetESDReader::~AliJetESDReader()
{
  // Destructor
    delete fChain;
    delete fTree;
    delete fESD;
    delete fTpcGrid;
    delete fEmcalGrid;
    if(fDZ)
      {
	delete fGrid0;
	delete fGrid1;
	delete fGrid2;
	delete fGrid3;
	delete fGrid4;
      }
}

//____________________________________________________________________________
void AliJetESDReader::OpenInputFiles()
{
  // chain for the ESDs
  fChain   = new TChain("esdTree");

  // get directory and pattern name from the header
   const char* dirName=fReaderHeader->GetDirectory();
   const char* pattern=fReaderHeader->GetPattern();
    
   // Add files matching patters to the chain
   void *dir  = gSystem->OpenDirectory(dirName);
   const char *name = 0x0;
   int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
   int a = 0;
   while ((name = gSystem->GetDirEntry(dir))){
       if (a>=nesd) continue;
       
       if (strstr(name,pattern)){
	 printf("Adding %s\n",name);
	 fChain->AddFile(Form("%s/%s/AliESDs.root",dirName,name));
	 a++;
       }
   }
  
  gSystem->FreeDirectory(dir);
  
  fESD = new AliESDEvent();
  fESD->ReadFromTree(fChain);

  int nMax = fChain->GetEntries(); 

  printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);

  // set number of events in header
  if (fReaderHeader->GetLastEvent() == -1)
    fReaderHeader->SetLastEvent(nMax);
  else {
    Int_t nUsr = fReaderHeader->GetLastEvent();
    fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
  }

}

//__________________________________________________________
void AliJetESDReader::SetApplyMIPCorrection(Bool_t val)
{
  //
  // Set flag to apply MIP correction fApplyMIPCorrection
  // - exclusive with fApplyFractionHadronicCorrection
  //

  fApplyMIPCorrection = val;
  if(fApplyMIPCorrection == kTRUE)
    {
      SetApplyFractionHadronicCorrection(kFALSE);
      printf("Enabling MIP Correction \n");
    }
  else 
    {
      printf("Disabling MIP Correction \n");
    }
}

//__________________________________________________________
void AliJetESDReader::SetApplyFractionHadronicCorrection(Bool_t val)
{
  //
  // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
  // - exclusive with fApplyMIPCorrection
  //

  fApplyFractionHadronicCorrection = val;
  if(fApplyFractionHadronicCorrection == kTRUE)
    {
      SetApplyMIPCorrection(kFALSE);
      printf("Enabling Fraction Hadronic Correction \n");
    }
  else 
    {
      printf("Disabling Fraction Hadronic Correction \n");
    }
}

//__________________________________________________________
void AliJetESDReader::SetFractionHadronicCorrection(Double_t val)
{
  //
  // Set value to fFractionHadronicCorrection (default is 0.3)
  // apply EMC hadronic correction fApplyFractionHadronicCorrection
  // - exclusive with fApplyMIPCorrection
  //

  fFractionHadronicCorrection = val;
  if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
    {
      SetApplyFractionHadronicCorrection(kTRUE);
      printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
    }
  else 
    {
      SetApplyFractionHadronicCorrection(kFALSE);
    }
}


//____________________________________________________________________________
void AliJetESDReader::SetInputEvent(const TObject* esd, const TObject* /*aod*/, const TObject* /*mc*/) {
    // Connect the tree
     fESD   = (AliESDEvent*) esd;
}

//____________________________________________________________________________
Bool_t AliJetESDReader::FillMomentumArray()
{
  // Fill momentum array

  Int_t goodTrack = 0;
  Int_t nt = 0;
  Float_t pt, eta;
  TVector3 p3;
  
  // clear momentum array
  ClearArray();
  fDebug = fReaderHeader->GetDebug();
  fOpt = fReaderHeader->GetDetector();

  if (!fESD) {
      return kFALSE;
  }

  // get number of tracks in event (for the loop)
  nt = fESD->GetNumberOfTracks();
  printf("Fill Momentum Array %5d  \n", nt);
 
  // temporary storage of signal and pt cut flag
  Int_t* sflag  = new Int_t[nt];
  Int_t* cflag  = new Int_t[nt];
  
  // get cuts set by user
  Float_t ptMin =  fReaderHeader->GetPtCut();
  Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
  Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
  
  //loop over tracks in ESD
  for (Int_t it = 0; it < nt; it++) {
      AliESDtrack *track = fESD->GetTrack(it);
      UInt_t status = track->GetStatus();
      
      Double_t mom[3];
      track->GetPxPyPz(mom);
      p3.SetXYZ(mom[0],mom[1],mom[2]);
      pt = p3.Pt();
      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();
      if ( (eta > etaMax) || (eta < etaMin))         continue;      // checking eta cut
      
      new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
      sflag[goodTrack]=0;
      if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
      cflag[goodTrack]=0;
      if (pt > ptMin) cflag[goodTrack]=1;                           // pt cut
      goodTrack++;
  }

  // set the signal flags
  fSignalFlag.Set(goodTrack,sflag);
  fCutFlag.Set(goodTrack,cflag);

  delete[] sflag;
  delete[] cflag;

  return kTRUE;

}

//____________________________________________________________________________
void AliJetESDReader::CreateTasks(TChain* tree)
{
  //
  // For reader task initialization
  //

  fDebug = fReaderHeader->GetDebug();
  fDZ = fReaderHeader->GetDZ();
  fTree = tree;

  // Init EMCAL geometry, if needed
   if(fGeom == 0)
	  SetEMCALGeometry();
	else Info(" SetEMCALGeometry","was already done.. it's called just once !!");

  // Init parameters
  InitParameters();
  // Create and init unit array
  InitUnitArray();

  // Create global reader task for analysis 
  fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
  // Create a task for to fill the charged particle information 
  fFillUAFromTracks = new AliJetESDFillUnitArrayTracks(); 
  fFillUAFromTracks->SetReaderHeader(fReaderHeader);
  fFillUAFromTracks->SetGeom(fGeom);
  fFillUAFromTracks->SetTPCGrid(fTpcGrid);
  fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
  if(fDZ)
    { // Calo dead zones inclusion
      fFillUAFromTracks->SetGrid0(fGrid0);
      fFillUAFromTracks->SetGrid1(fGrid1);
      fFillUAFromTracks->SetGrid2(fGrid2);
      fFillUAFromTracks->SetGrid3(fGrid3);
      fFillUAFromTracks->SetGrid4(fGrid4);
    }
  fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
  fFillUAFromTracks->SetHadCorrector(fHadCorr);
  // Create a task for to fill the neutral particle information 
  fFillUAFromEMCalDigits = new AliJetESDFillUnitArrayEMCalDigits();
  fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
  fFillUAFromEMCalDigits->SetGeom(fGeom);
  fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
  fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
  fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
  fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
   fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
  // Add the task to global task
  fFillUnitArray->Add(fFillUAFromTracks);
  fFillUnitArray->Add(fFillUAFromEMCalDigits);
  fFillUAFromTracks->SetActive(kFALSE);
  fFillUAFromEMCalDigits->SetActive(kFALSE);

  cout << "Tasks instantiated at that stage ! " << endl;
  cout << "You can loop over events now ! " << endl;
   
}

//____________________________________________________________________________
Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
{
  //
  // Reader task execussion
  //

  fProcId = procid;
  fRefArray = refArray;

  // clear momentum array
  ClearArray();

  fDebug = fReaderHeader->GetDebug();
  fOpt = fReaderHeader->GetDetector();

  if(!fESD) {
    return kFALSE;
  }
  
  // TPC only or Digits+TPC or Clusters+TPC
  if(fOpt%2==!0 && fOpt!=0){ 
    fFillUAFromTracks->SetESD(fESD);
    fFillUAFromTracks->SetActive(kTRUE);
    fFillUAFromTracks->SetUnitArray(fUnitArray);
    fFillUAFromTracks->SetRefArray(fRefArray);
    fFillUAFromTracks->SetProcId(fProcId);
    //    fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed
    fFillUAFromTracks->Exec("tpc");              // Temporarily added
    if(fOpt==1){
      fNumCandidate = fFillUAFromTracks->GetMult();
      fNumCandidateCut = fFillUAFromTracks->GetMultCut();
    }
  }

  // Digits only or Digits+TPC  
  if(fOpt>=2 && fOpt<=3){
    fFillUAFromEMCalDigits->SetESD(fESD);
    fFillUAFromEMCalDigits->SetActive(kTRUE);
    fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
    fFillUAFromEMCalDigits->SetRefArray(fRefArray);
    fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
    fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
    fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
    fFillUAFromEMCalDigits->Exec("digits");      // Temporarily added
    fNumCandidate = fFillUAFromEMCalDigits->GetMult();
    fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
  }

  //  fFillUnitArray->ExecuteTask();             // Temporarily commented

  return kTRUE;
}


//____________________________________________________________________________  
void AliJetESDReader::InitParameters()
{
  // Initialise parameters
  fOpt = fReaderHeader->GetDetector();
  //  fHadCorr        = 0;          // For hadron correction
  if(fEFlag==kFALSE){
    if(fOpt==0 || fOpt==1)  
      fECorrection    = 0;        // For electron correction
    else fECorrection = 1;        // For electron correction
  }
  fNumUnits       = fGeom->GetEMCGeometry()->GetNCells();      // Number of cells in EMCAL
  if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
}

//____________________________________________________________________________
void AliJetESDReader::InitUnitArray()
{
  //
  // Create and Initialise unit arrays
  //

  Int_t nElements = fTpcGrid->GetNEntries();
  Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
  if(fArrayInitialised) fUnitArray->Delete();

  if(fTpcGrid->GetGridType()==0)
    { // Fill the following quantities :
      // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi, 
      // detector flag, in/out jet, pt cut, mass, cluster ID)
      for(Int_t nBin = 1; nBin < nElements+1; nBin++)
	{
	  fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
	  phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
	  Deta = fTpcGrid->GetDeta();
	  Dphi = fTpcGrid->GetDphi();
	  new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	}
    }

  if(fTpcGrid->GetGridType()==1)
    {
      Int_t nGaps = 0;
      Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;

      if(fDZ)
	{
	  // Define a grid of cell for the gaps between SM
	  Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
	  Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
	  fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
	  fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
	  fGrid0->SetGridType(0);
	  fGrid0->SetMatrixIndexes();
	  fGrid0->SetIndexIJ();
	  n0 = fGrid0->GetNEntries();
	  fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
	  fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
	  fGrid1->SetGridType(0);
	  fGrid1->SetMatrixIndexes();
	  fGrid1->SetIndexIJ();
	  n1 = fGrid1->GetNEntries();
	  fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
	  fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
	  fGrid2->SetGridType(0);
	  fGrid2->SetMatrixIndexes();
	  fGrid2->SetIndexIJ();
	  n2 = fGrid2->GetNEntries();
	  fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
	  fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
	  fGrid3->SetGridType(0);  
	  fGrid3->SetMatrixIndexes();
	  fGrid3->SetIndexIJ();
	  n3 = fGrid3->GetNEntries();
	  fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
	  fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
	  fGrid4->SetGridType(0);
	  fGrid4->SetMatrixIndexes();
	  fGrid4->SetIndexIJ();
	  n4 = fGrid4->GetNEntries();

	  nGaps = n0+n1+n2+n3+n4;

	}

      for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++) 
	{
	  if(nBin<fNumUnits)
	    {
	      fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry 
	      // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
	      phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
	      Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
	      Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
	      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	    } 
	  else {
	    if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
	      fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
	      phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
	      Deta = fTpcGrid->GetDeta();
	      Dphi = fTpcGrid->GetDphi();
	      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	    }
	    else {
	      if(fDZ) {
		if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
		  if(nBin<fNumUnits+nElements+n0)
		    {
		      Float_t phi0 = eta = 0.;
		      fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi0,eta);
		      Deta = fGrid0->GetDeta(); 
		      Dphi = fGrid0->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		  else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
		    {
		      Float_t phi0 = eta = 0.;
		      fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi0,eta);
		      Deta = fGrid1->GetDeta(); 
		      Dphi = fGrid1->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		  else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
		    {
		      Float_t phi0 = eta = 0.;
		      fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi0,eta);
		      Deta = fGrid2->GetDeta(); 
		      Dphi = fGrid2->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		  else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
		    {
		      Float_t phi0 = eta = 0.;
		      fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi0,eta);
		      Deta = fGrid3->GetDeta(); 
		      Dphi = fGrid3->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		  else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
		    {
		      Float_t phi0 = eta = 0.;
		      fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi0,eta);
		      Deta = fGrid4->GetDeta(); 
		      Dphi = fGrid4->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		}
	      } // end if(fDZ)
	    } // end else 2
	  } // end else 1
	} // end loop on nBin
    } // end grid type == 1
  fArrayInitialised = 1;
}
Commit	Line	Data
	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	// Jet ESD Reader
	18	// ESD reader for jet analysis
	19	// Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
	20	// Magali Estienne <magali.estienne@subatech.in2p3.fr>
	21	//
	22	// **February 2011
	23	// implemented standard geometry (AliEMCALGeometry) instead of dummy one (AliJetDummyGeo)
	24	// moved geometry definition in AliJetReader
	25	// marco.bregant@subatech.in2p3.fr
	26	//-------------------------------------------------------------------------
	27
	28	// --- Standard library ---
	29	#include <Riostream.h>
	30
	31	// --- ROOT system ---
	32	#include <TSystem.h>
	33	#include <TString.h>
	34	#include <TStopwatch.h>
	35	#include <TLorentzVector.h>
	36	#include <TVector3.h>
	37	#include "TTask.h"
	38	#include "TTree.h"
	39	#include "TFile.h"
	40	#include <TGeoManager.h>
	41	#include <assert.h>
	42	#include <TRefArray.h>
	43	#include <TMath.h>
	44	#include <TProcessID.h>
	45	#include <TRandom3.h>
	46
	47	// --- AliRoot header files ---
	48	#include "AliJetESDReader.h"
	49	#include "AliJetESDReaderHeader.h"
	50	#include "AliESDEvent.h"
	51	#include "AliVEvent.h"
	52	#include "AliESD.h"
	53	#include "AliESDtrack.h"
	54	#include "AliEMCALGeometry.h"
	55	#include "AliEMCALEMCGeometry.h"
	56	#include "AliJetESDFillUnitArrayTracks.h"
	57	#include "AliJetESDFillUnitArrayEMCalDigits.h"
	58	#include "AliJetHadronCorrectionv1.h"
	59	#include "AliJetUnitArray.h"
	60	#include "AliAnalysisTask.h"
	61	#include "AliOADBContainer.h"
	62
	63	using std::cout;
	64	using std::endl;
	65	ClassImp(AliJetESDReader)
	66
	67	AliJetESDReader::AliJetESDReader():
	68	AliJetReader(),
	69	fHadCorr(0x0),
	70	fTpcGrid(0x0),
	71	fEmcalGrid(0x0),
	72	fGrid0(0),
	73	fGrid1(0),
	74	fGrid2(0),
	75	fGrid3(0),
	76	fGrid4(0),
	77	fApplyElectronCorrection(kFALSE),
	78	fApplyMIPCorrection(kTRUE),
	79	fApplyFractionHadronicCorrection(kFALSE),
	80	fFractionHadronicCorrection(0.3),
	81	fNumUnits(0),
	82	fDebug(0),
	83	fMass(0),
	84	fSign(0),
	85	fNIn(0),
	86	fOpt(0),
	87	fDZ(0),
	88	fNeta(0),
	89	fNphi(0),
	90	fRefArray(0x0),
	91	fProcId(kFALSE)
	92	{
	93	// Constructor
	94	}
	95
	96	//____________________________________________________________________________
	97	AliJetESDReader::~AliJetESDReader()
	98	{
	99	// Destructor
	100	delete fChain;
	101	delete fTree;
	102	delete fESD;
	103	delete fTpcGrid;
	104	delete fEmcalGrid;
	105	if(fDZ)
	106	{
	107	delete fGrid0;
	108	delete fGrid1;
	109	delete fGrid2;
	110	delete fGrid3;
	111	delete fGrid4;
	112	}
	113	}
	114
	115	//____________________________________________________________________________
	116	void AliJetESDReader::OpenInputFiles()
	117	{
	118	// chain for the ESDs
	119	fChain = new TChain("esdTree");
	120
	121	// get directory and pattern name from the header
	122	const char* dirName=fReaderHeader->GetDirectory();
	123	const char* pattern=fReaderHeader->GetPattern();
	124
	125	// Add files matching patters to the chain
	126	void *dir = gSystem->OpenDirectory(dirName);
	127	const char *name = 0x0;
	128	int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
	129	int a = 0;
	130	while ((name = gSystem->GetDirEntry(dir))){
	131	if (a>=nesd) continue;
	132
	133	if (strstr(name,pattern)){
	134	printf("Adding %s\n",name);
	135	fChain->AddFile(Form("%s/%s/AliESDs.root",dirName,name));
	136	a++;
	137	}
	138	}
	139
	140	gSystem->FreeDirectory(dir);
	141
	142	fESD = new AliESDEvent();
	143	fESD->ReadFromTree(fChain);
	144
	145	int nMax = fChain->GetEntries();
	146
	147	printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
	148
	149	// set number of events in header
	150	if (fReaderHeader->GetLastEvent() == -1)
	151	fReaderHeader->SetLastEvent(nMax);
	152	else {
	153	Int_t nUsr = fReaderHeader->GetLastEvent();
	154	fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
	155	}
	156
	157	}
	158
	159	//__________________________________________________________
	160	void AliJetESDReader::SetApplyMIPCorrection(Bool_t val)
	161	{
	162	//
	163	// Set flag to apply MIP correction fApplyMIPCorrection
	164	// - exclusive with fApplyFractionHadronicCorrection
	165	//
	166
	167	fApplyMIPCorrection = val;
	168	if(fApplyMIPCorrection == kTRUE)
	169	{
	170	SetApplyFractionHadronicCorrection(kFALSE);
	171	printf("Enabling MIP Correction \n");
	172	}
	173	else
	174	{
	175	printf("Disabling MIP Correction \n");
	176	}
	177	}
	178
	179	//__________________________________________________________
	180	void AliJetESDReader::SetApplyFractionHadronicCorrection(Bool_t val)
	181	{
	182	//
	183	// Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
	184	// - exclusive with fApplyMIPCorrection
	185	//
	186
	187	fApplyFractionHadronicCorrection = val;
	188	if(fApplyFractionHadronicCorrection == kTRUE)
	189	{
	190	SetApplyMIPCorrection(kFALSE);
	191	printf("Enabling Fraction Hadronic Correction \n");
	192	}
	193	else
	194	{
	195	printf("Disabling Fraction Hadronic Correction \n");
	196	}
	197	}
	198
	199	//__________________________________________________________
	200	void AliJetESDReader::SetFractionHadronicCorrection(Double_t val)
	201	{
	202	//
	203	// Set value to fFractionHadronicCorrection (default is 0.3)
	204	// apply EMC hadronic correction fApplyFractionHadronicCorrection
	205	// - exclusive with fApplyMIPCorrection
	206	//
	207
	208	fFractionHadronicCorrection = val;
	209	if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
	210	{
	211	SetApplyFractionHadronicCorrection(kTRUE);
	212	printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
	213	}
	214	else
	215	{
	216	SetApplyFractionHadronicCorrection(kFALSE);
	217	}
	218	}
	219
	220
	221	//____________________________________________________________________________
	222	void AliJetESDReader::SetInputEvent(const TObject* esd, const TObject* /aod/, const TObject* /mc/) {
	223	// Connect the tree
	224	fESD = (AliESDEvent*) esd;
	225	}
	226
	227	//____________________________________________________________________________
	228	Bool_t AliJetESDReader::FillMomentumArray()
	229	{
	230	// Fill momentum array
	231
	232	Int_t goodTrack = 0;
	233	Int_t nt = 0;
	234	Float_t pt, eta;
	235	TVector3 p3;
	236
	237	// clear momentum array
	238	ClearArray();
	239	fDebug = fReaderHeader->GetDebug();
	240	fOpt = fReaderHeader->GetDetector();
	241
	242	if (!fESD) {
	243	return kFALSE;
	244	}
	245
	246	// get number of tracks in event (for the loop)
	247	nt = fESD->GetNumberOfTracks();
	248	printf("Fill Momentum Array %5d \n", nt);
	249
	250	// temporary storage of signal and pt cut flag
	251	Int_t* sflag = new Int_t[nt];
	252	Int_t* cflag = new Int_t[nt];
	253
	254	// get cuts set by user
	255	Float_t ptMin = fReaderHeader->GetPtCut();
	256	Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
	257	Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
	258
	259	//loop over tracks in ESD
	260	for (Int_t it = 0; it < nt; it++) {
	261	AliESDtrack *track = fESD->GetTrack(it);
	262	UInt_t status = track->GetStatus();
	263
	264	Double_t mom[3];
	265	track->GetPxPyPz(mom);
	266	p3.SetXYZ(mom[0],mom[1],mom[2]);
	267	pt = p3.Pt();
	268	if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
	269	if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
	270	if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
	271	&& TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
	272	if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
	273	&& TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
	274	eta = p3.Eta();
	275	if ( (eta > etaMax) \|\| (eta < etaMin)) continue; // checking eta cut
	276
	277	new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
	278	sflag[goodTrack]=0;
	279	if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
	280	cflag[goodTrack]=0;
	281	if (pt > ptMin) cflag[goodTrack]=1; // pt cut
	282	goodTrack++;
	283	}
	284
	285	// set the signal flags
	286	fSignalFlag.Set(goodTrack,sflag);
	287	fCutFlag.Set(goodTrack,cflag);
	288
	289	delete[] sflag;
	290	delete[] cflag;
	291
	292	return kTRUE;
	293
	294	}
	295
	296	//____________________________________________________________________________
	297	void AliJetESDReader::CreateTasks(TChain* tree)
	298	{
	299	//
	300	// For reader task initialization
	301	//
	302
	303	fDebug = fReaderHeader->GetDebug();
	304	fDZ = fReaderHeader->GetDZ();
	305	fTree = tree;
	306
	307	// Init EMCAL geometry, if needed
	308	if(fGeom == 0)
	309	SetEMCALGeometry();
	310	else Info(" SetEMCALGeometry","was already done.. it's called just once !!");
	311
	312	// Init parameters
	313	InitParameters();
	314	// Create and init unit array
	315	InitUnitArray();
	316
	317	// Create global reader task for analysis
	318	fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
	319	// Create a task for to fill the charged particle information
	320	fFillUAFromTracks = new AliJetESDFillUnitArrayTracks();
	321	fFillUAFromTracks->SetReaderHeader(fReaderHeader);
	322	fFillUAFromTracks->SetGeom(fGeom);
	323	fFillUAFromTracks->SetTPCGrid(fTpcGrid);
	324	fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
	325	if(fDZ)
	326	{ // Calo dead zones inclusion
	327	fFillUAFromTracks->SetGrid0(fGrid0);
	328	fFillUAFromTracks->SetGrid1(fGrid1);
	329	fFillUAFromTracks->SetGrid2(fGrid2);
	330	fFillUAFromTracks->SetGrid3(fGrid3);
	331	fFillUAFromTracks->SetGrid4(fGrid4);
	332	}
	333	fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
	334	fFillUAFromTracks->SetHadCorrector(fHadCorr);
	335	// Create a task for to fill the neutral particle information
	336	fFillUAFromEMCalDigits = new AliJetESDFillUnitArrayEMCalDigits();
	337	fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
	338	fFillUAFromEMCalDigits->SetGeom(fGeom);
	339	fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
	340	fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
	341	fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
	342	fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
	343	fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
	344	// Add the task to global task
	345	fFillUnitArray->Add(fFillUAFromTracks);
	346	fFillUnitArray->Add(fFillUAFromEMCalDigits);
	347	fFillUAFromTracks->SetActive(kFALSE);
	348	fFillUAFromEMCalDigits->SetActive(kFALSE);
	349
	350	cout << "Tasks instantiated at that stage ! " << endl;
	351	cout << "You can loop over events now ! " << endl;
	352
	353	}
	354
	355	//____________________________________________________________________________
	356	Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
	357	{
	358	//
	359	// Reader task execussion
	360	//
	361
	362	fProcId = procid;
	363	fRefArray = refArray;
	364
	365	// clear momentum array
	366	ClearArray();
	367
	368	fDebug = fReaderHeader->GetDebug();
	369	fOpt = fReaderHeader->GetDetector();
	370
	371	if(!fESD) {
	372	return kFALSE;
	373	}
	374
	375	// TPC only or Digits+TPC or Clusters+TPC
	376	if(fOpt%2==!0 && fOpt!=0){
	377	fFillUAFromTracks->SetESD(fESD);
	378	fFillUAFromTracks->SetActive(kTRUE);
	379	fFillUAFromTracks->SetUnitArray(fUnitArray);
	380	fFillUAFromTracks->SetRefArray(fRefArray);
	381	fFillUAFromTracks->SetProcId(fProcId);
	382	// fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed
	383	fFillUAFromTracks->Exec("tpc"); // Temporarily added
	384	if(fOpt==1){
	385	fNumCandidate = fFillUAFromTracks->GetMult();
	386	fNumCandidateCut = fFillUAFromTracks->GetMultCut();
	387	}
	388	}
	389
	390	// Digits only or Digits+TPC
	391	if(fOpt>=2 && fOpt<=3){
	392	fFillUAFromEMCalDigits->SetESD(fESD);
	393	fFillUAFromEMCalDigits->SetActive(kTRUE);
	394	fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
	395	fFillUAFromEMCalDigits->SetRefArray(fRefArray);
	396	fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
	397	fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
	398	fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
	399	fFillUAFromEMCalDigits->Exec("digits"); // Temporarily added
	400	fNumCandidate = fFillUAFromEMCalDigits->GetMult();
	401	fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
	402	}
	403
	404	// fFillUnitArray->ExecuteTask(); // Temporarily commented
	405
	406	return kTRUE;
	407	}
	408
	409
	410	//____________________________________________________________________________
	411	void AliJetESDReader::InitParameters()
	412	{
	413	// Initialise parameters
	414	fOpt = fReaderHeader->GetDetector();
	415	// fHadCorr = 0; // For hadron correction
	416	if(fEFlag==kFALSE){
	417	if(fOpt==0 \|\| fOpt==1)
	418	fECorrection = 0; // For electron correction
	419	else fECorrection = 1; // For electron correction
	420	}
	421	fNumUnits = fGeom->GetEMCGeometry()->GetNCells(); // Number of cells in EMCAL
	422	if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
	423	}
	424
	425	//____________________________________________________________________________
	426	void AliJetESDReader::InitUnitArray()
	427	{
	428	//
	429	// Create and Initialise unit arrays
	430	//
	431
	432	Int_t nElements = fTpcGrid->GetNEntries();
	433	Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
	434	if(fArrayInitialised) fUnitArray->Delete();
	435
	436	if(fTpcGrid->GetGridType()==0)
	437	{ // Fill the following quantities :
	438	// Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi,
	439	// detector flag, in/out jet, pt cut, mass, cluster ID)
	440	for(Int_t nBin = 1; nBin < nElements+1; nBin++)
	441	{
	442	fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
	443	phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
	444	Deta = fTpcGrid->GetDeta();
	445	Dphi = fTpcGrid->GetDphi();
	446	new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	447	}
	448	}
	449
	450	if(fTpcGrid->GetGridType()==1)
	451	{
	452	Int_t nGaps = 0;
	453	Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
	454
	455	if(fDZ)
	456	{
	457	// Define a grid of cell for the gaps between SM
	458	Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
	459	Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
	460	fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
	461	fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
	462	fGrid0->SetGridType(0);
	463	fGrid0->SetMatrixIndexes();
	464	fGrid0->SetIndexIJ();
	465	n0 = fGrid0->GetNEntries();
	466	fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
	467	fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
	468	fGrid1->SetGridType(0);
	469	fGrid1->SetMatrixIndexes();
	470	fGrid1->SetIndexIJ();
	471	n1 = fGrid1->GetNEntries();
	472	fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
	473	fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
	474	fGrid2->SetGridType(0);
	475	fGrid2->SetMatrixIndexes();
	476	fGrid2->SetIndexIJ();
	477	n2 = fGrid2->GetNEntries();
	478	fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
	479	fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
	480	fGrid3->SetGridType(0);
	481	fGrid3->SetMatrixIndexes();
	482	fGrid3->SetIndexIJ();
	483	n3 = fGrid3->GetNEntries();
	484	fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
	485	fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
	486	fGrid4->SetGridType(0);
	487	fGrid4->SetMatrixIndexes();
	488	fGrid4->SetIndexIJ();
	489	n4 = fGrid4->GetNEntries();
	490
	491	nGaps = n0+n1+n2+n3+n4;
	492
	493	}
	494
	495	for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
	496	{
	497	if(nBin<fNumUnits)
	498	{
	499	fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
	500	// fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
	501	phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
	502	Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
	503	Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
	504	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	505	}
	506	else {
	507	if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
	508	fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
	509	phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
	510	Deta = fTpcGrid->GetDeta();
	511	Dphi = fTpcGrid->GetDphi();
	512	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	513	}
	514	else {
	515	if(fDZ) {
	516	if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
	517	if(nBin<fNumUnits+nElements+n0)
	518	{
	519	Float_t phi0 = eta = 0.;
	520	fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi0,eta);
	521	Deta = fGrid0->GetDeta();
	522	Dphi = fGrid0->GetDphi();
	523	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	524	}
	525	else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
	526	{
	527	Float_t phi0 = eta = 0.;
	528	fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi0,eta);
	529	Deta = fGrid1->GetDeta();
	530	Dphi = fGrid1->GetDphi();
	531	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	532	}
	533	else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
	534	{
	535	Float_t phi0 = eta = 0.;
	536	fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi0,eta);
	537	Deta = fGrid2->GetDeta();
	538	Dphi = fGrid2->GetDphi();
	539	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	540	}
	541	else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
	542	{
	543	Float_t phi0 = eta = 0.;
	544	fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi0,eta);
	545	Deta = fGrid3->GetDeta();
	546	Dphi = fGrid3->GetDphi();
	547	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	548	}
	549	else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
	550	{
	551	Float_t phi0 = eta = 0.;
	552	fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi0,eta);
	553	Deta = fGrid4->GetDeta();
	554	Dphi = fGrid4->GetDphi();
	555	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	556	}
	557	}
	558	} // end if(fDZ)
	559	} // end else 2
	560	} // end else 1
	561	} // end loop on nBin
	562	} // end grid type == 1
	563	fArrayInitialised = 1;
	564	}
	565
	566
	567