[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>
//------------------------------------------------------------------------- 

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

// --- ROOT system ---
#include <TSystem.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 "AliESD.h"
#include "AliESDtrack.h"
#include "AliJetDummyGeo.h"
#include "AliJetESDFillUnitArrayTracks.h"
#include "AliJetESDFillUnitArrayEMCalDigits.h"
#include "AliJetUnitArray.h"
#include "AliAnalysisTask.h"

ClassImp(AliJetESDReader)

AliJetESDReader::AliJetESDReader():
  AliJetReader(),  
  fGeom(0),
  fHadCorr(0x0),
  fTpcGrid(0x0),
  fEmcalGrid(0x0),
  fGrid0(0),
  fGrid1(0),
  fGrid2(0),
  fGrid3(0),
  fGrid4(0),
  fPtCut(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);
	 char path[256];
	 sprintf(path,"%s/%s/AliESDs.root",dirName,name);
	 fChain->AddFile(path);
	 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 
  SetEMCALGeometry();
  // 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;
}

//____________________________________________________________________________
Bool_t AliJetESDReader::SetEMCALGeometry()
{
  // 
  // Set the EMCal Geometry
  //

  if (!fTree->GetFile()) 
    return kFALSE;

  TString geomFile(fTree->GetFile()->GetName());
  geomFile.ReplaceAll("AliESDs", "geometry");
  
  // temporary workaround for PROOF bug #18505
  geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
  if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());

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