[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 "AliJetFillUnitArrayTracks.h"
#include "AliJetFillUnitArrayEMCalDigits.h"
#include "AliJetUnitArray.h"
#include "AliAnalysisTask.h"

ClassImp(AliJetESDReader)

AliJetESDReader::AliJetESDReader():
  AliJetReader(),  
  fGeom(0),
  fChain(0x0),
  fTree(0x0),
  fESD(0x0),
  fHadCorr(0x0),
  fTpcGrid(0x0),
  fEmcalGrid(0x0),
  fGrid0(0),
  fGrid1(0),
  fGrid2(0),
  fGrid3(0),
  fGrid4(0),
  fPtCut(0),
  fHCorrection(0),
  fECorrection(0),
  fEFlag(kFALSE),
  fNumUnits(0),
  fDebug(0),
  fMass(0),
  fSign(0),
  fNIn(0),
  fOpt(0),
  fDZ(0),
  fNeta(0),
  fNphi(0),
  fArrayInitialised(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::SetInputEvent(TObject* esd, TObject* /*aod*/, 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 AliJetFillUnitArrayTracks(); 
  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->SetHadCorrection(fHCorrection);
  fFillUAFromTracks->SetHadCorrector(fHadCorr);
  // Create a task for to fill the neutral particle information 
  fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits();
  fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
  fFillUAFromEMCalDigits->SetGeom(fGeom);
  fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
  fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
  fFillUAFromEMCalDigits->SetEleCorrection(fECorrection);
  // 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;
  vector<Float_t> vtmp(3);

  // 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();

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