[u/mrichter/AliRoot.git] / JETAN / AliJetAODReader.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 AOD Reader 
// AOD reader for jet analysis
// This is the reader which must be used if the jet analysis task
// is executed after the ESD filter task, in order to read its output
//
// Author: Davide Perrino <davide.perrino@cern.ch>
//------------------------------------------------------------------------- 


#include <Riostream.h>
#include <TSystem.h>
#include <TLorentzVector.h>
#include <TVector3.h>
#include <TChain.h>
#include <TFile.h>
#include <TTask.h>
#include <TGeoManager.h>

#include "AliJetAODReader.h"
#include "AliJetAODReaderHeader.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliJetDummyGeo.h"
#include "AliJetAODFillUnitArrayTracks.h"
#include "AliJetAODFillUnitArrayEMCalDigits.h"
#include "AliJetHadronCorrection.h"
#include "AliJetUnitArray.h"

ClassImp(AliJetAODReader)

AliJetAODReader::AliJetAODReader():
    AliJetReader(),
    fChain(0x0),
    fTree(0x0),
    fAOD(0x0),
    fRef(new TRefArray),
    fDebug(0),
    fOpt(0),
    fGeom(0),
    fHadCorr(0x0),
    fTpcGrid(0x0),
    fEmcalGrid(0x0),
    fGrid0(0),
    fGrid1(0),
    fGrid2(0),
    fGrid3(0),
    fGrid4(0),
    fPtCut(0),
    fHCorrection(0),
    fApplyElectronCorrection(kFALSE),
    fEFlag(kFALSE),
    fApplyMIPCorrection(kTRUE),
    fApplyFractionHadronicCorrection(kFALSE),
    fFractionHadronicCorrection(0.3),
    fNumUnits(0),
    fMass(0),
    fSign(0),
    fNIn(0),
    fDZ(0),
    fNeta(0),
    fNphi(0),
    fArrayInitialised(0),
    fRefArray(0x0),
    fProcId(kFALSE)

{
  // Constructor    
}

//____________________________________________________________________________

AliJetAODReader::~AliJetAODReader()
{
  // Destructor
    delete fChain;
    delete fAOD;
    delete fRef;
    delete fTpcGrid;
    delete fEmcalGrid;
    if(fDZ)
      {
        delete fGrid0;
        delete fGrid1;
        delete fGrid2;
        delete fGrid3;
        delete fGrid4;
      }

}

//____________________________________________________________________________

void AliJetAODReader::OpenInputFiles()
{
    // Open the necessary input files
    // chain for the AODs
  fChain   = new TChain("aodTree");

  // 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 naod = ((AliJetAODReaderHeader*) fReaderHeader)->GetNaod();
   int a = 0;
   while ((name = gSystem->GetDirEntry(dir))){
       if (a>=naod) continue;
       
       if (strstr(name,pattern)){
	   char path[256];
	   sprintf(path,"%s/%s/aod.root",dirName,name);
	   fChain->AddFile(path);
	   a++;
       }
   }
  
  gSystem->FreeDirectory(dir);
  
  fAOD = 0;
  fChain->SetBranchAddress("AOD",&fAOD);
  
  int nMax = fChain->GetEntries(); 

  printf("\n AliJetAODReader: 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 AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) {
    // Connect the tree
    // For AOD reader it's needed only to set the number of events
     fChain = (TChain*)      tree;
     
     Int_t nMax = fChain->GetEntries(); 
     printf("\n AliJetAODReader: Total number of events in chain= %5d \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));
     }
}

//____________________________________________________________________________

Bool_t AliJetAODReader::FillMomentumArray()
{
  // Clear momentum array
  ClearArray();
  fRef->Clear();
  fDebug = fReaderHeader->GetDebug();
  
  if (!fAOD) {
      return kFALSE;
  }
  
  // get number of tracks in event (for the loop)
  Int_t nt = fAOD->GetNTracks();
  printf("AOD tracks: %5d \t", 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();  
  UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();

  //loop over tracks
  Int_t aodTrack = 0;
  Float_t pt, eta;
  TVector3 p3;

  for (Int_t it = 0; it < nt; it++) {
    AliAODTrack *track = fAOD->GetTrack(it);
    UInt_t status = track->GetStatus();
    
    Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
    p3.SetXYZ(mom[0],mom[1],mom[2]);
    pt = p3.Pt();
    eta = p3.Eta();
    if (status == 0) continue;
    if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
    if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut

    new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
    sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
    cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
    aodTrack++;
    fRef->Add(track);
  }
  printf("Used AOD tracks: %5d \n", aodTrack);
  // set the signal flags
  fSignalFlag.Set(aodTrack,sflag);
  fCutFlag.Set(aodTrack,cflag);

  delete [] sflag;
  delete [] cflag;
  
  return kTRUE;
}

//__________________________________________________________
void AliJetAODReader::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 AliJetAODReader::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 AliJetAODReader::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 AliJetAODReader::CreateTasks(TChain* tree)
{
  fDebug = fReaderHeader->GetDebug();
  fDZ = fReaderHeader->GetDZ();
  fTree = tree;

  // Init EMCAL geometry and create UnitArray object
  SetEMCALGeometry();
  //  cout << "In create task" << endl;
  InitParameters();
  InitUnitArray();

  fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
  fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
  fFillUAFromTracks->SetReaderHeader(fReaderHeader);
  fFillUAFromTracks->SetGeom(fGeom);
  fFillUAFromTracks->SetTPCGrid(fTpcGrid);
  fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);

  if(fDZ)
    {
      fFillUAFromTracks->SetGrid0(fGrid0);
      fFillUAFromTracks->SetGrid1(fGrid1);
      fFillUAFromTracks->SetGrid2(fGrid2);
      fFillUAFromTracks->SetGrid3(fGrid3);
      fFillUAFromTracks->SetGrid4(fGrid4);
    }
  fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
  fFillUAFromTracks->SetHadCorrector(fHadCorr);
  fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
  fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
  fFillUAFromEMCalDigits->SetGeom(fGeom);
  fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
  fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
  fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
  fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
  fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);

  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 AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
{
  //
  // Main function
  // Fill the reader part
  //

  fProcId = procid;
  fRefArray = refArray;
//(not used ?)  Int_t nEntRef = fRefArray->GetEntries();
//(not used ?)  Int_t nEntUnit = fUnitArray->GetEntries();

  // clear momentum array
  ClearArray();

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

  if(!fAOD) {
    return kFALSE;
  }

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

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