[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(),
    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),
    fApplyElectronCorrection(kFALSE),
    fApplyMIPCorrection(kTRUE),
    fApplyFractionHadronicCorrection(kFALSE),
    fFractionHadronicCorrection(0.3),
    fNumUnits(0),
    fMass(0),
    fSign(0),
    fNIn(0),
    fDZ(0),
    fNeta(0),
    fNphi(0),
    fRefArray(0x0),
    fProcId(kFALSE)
{
  // Constructor    
}

//____________________________________________________________________________

AliJetAODReader::~AliJetAODReader()
{
  // Destructor
    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

    if (aodTrack == 0) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
    
    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)
{
  //
  // For reader task initialization
  //

  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., deltaEta = 0., deltaPhi = 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, deltaEta, deltaPhi,
      // 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);
          deltaEta = fTpcGrid->GetDeta();
          deltaPhi = fTpcGrid->GetDphi();
          new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,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);
	      deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
	      deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
	      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,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);
	      deltaEta = fTpcGrid->GetDeta();
	      deltaPhi = fTpcGrid->GetDphi();
	      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	    }
	    else {
	      if(fDZ) {
		if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
		  if(nBin<fNumUnits+nElements+n0)
		    {
		      phi = eta = 0.;
		      fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
		      deltaEta = fGrid0->GetDeta(); 
		      deltaPhi = fGrid0->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		  else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
		    {
		      phi = eta = 0.;
		      fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
		      deltaEta = fGrid1->GetDeta(); 
		      deltaPhi = fGrid1->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		  else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
		    {
		      phi = eta = 0.;
		      fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
		      deltaEta = fGrid2->GetDeta(); 
		      deltaPhi = fGrid2->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		  else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
		    {
		      phi = eta = 0.;
		      fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
		      deltaEta = fGrid3->GetDeta(); 
		      deltaPhi = fGrid3->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
		    }
		  else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
		    {
		      phi = eta = 0.;
		      fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
		      deltaEta = fGrid4->GetDeta(); 
		      deltaPhi = fGrid4->GetDphi(); 
		      new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,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(),
586f2bc3	49	fAOD(0x0),
	50	fRef(new TRefArray),
	51	fDebug(0),
be6e5811	52	fOpt(0),
	53	fGeom(0),
	54	fHadCorr(0x0),
	55	fTpcGrid(0x0),
	56	fEmcalGrid(0x0),
	57	fGrid0(0),
	58	fGrid1(0),
	59	fGrid2(0),
	60	fGrid3(0),
	61	fGrid4(0),
	62	fPtCut(0),
be6e5811	63	fApplyElectronCorrection(kFALSE),
be6e5811	64	fApplyMIPCorrection(kTRUE),
	65	fApplyFractionHadronicCorrection(kFALSE),
	66	fFractionHadronicCorrection(0.3),
	67	fNumUnits(0),
	68	fMass(0),
	69	fSign(0),
	70	fNIn(0),
	71	fDZ(0),
	72	fNeta(0),
	73	fNphi(0),
be6e5811	74	fRefArray(0x0),
be6e5811	75	fProcId(kFALSE)
586f2bc3	76	{
	77	// Constructor
	78	}
	79
	80	//____________________________________________________________________________
	81
	82	AliJetAODReader::~AliJetAODReader()
	83	{
	84	// Destructor
586f2bc3	85	delete fAOD;
586f2bc3	86	delete fRef;
be6e5811	87	delete fTpcGrid;
	88	delete fEmcalGrid;
	89	if(fDZ)
	90	{
	91	delete fGrid0;
	92	delete fGrid1;
	93	delete fGrid2;
	94	delete fGrid3;
	95	delete fGrid4;
	96	}
	97
586f2bc3	98	}
	99
	100	//____________________________________________________________________________
	101
	102	void AliJetAODReader::OpenInputFiles()
	103	{
	104	// Open the necessary input files
	105	// chain for the AODs
	106	fChain = new TChain("aodTree");
	107
	108	// get directory and pattern name from the header
	109	const char* dirName=fReaderHeader->GetDirectory();
	110	const char* pattern=fReaderHeader->GetPattern();
	111
	112	// // Add files matching patters to the chain
	113
	114	void *dir = gSystem->OpenDirectory(dirName);
	115	const char *name = 0x0;
	116	int naod = ((AliJetAODReaderHeader*) fReaderHeader)->GetNaod();
	117	int a = 0;
	118	while ((name = gSystem->GetDirEntry(dir))){
	119	if (a>=naod) continue;
	120
	121	if (strstr(name,pattern)){
	122	char path[256];
	123	sprintf(path,"%s/%s/aod.root",dirName,name);
	124	fChain->AddFile(path);
	125	a++;
	126	}
	127	}
	128
	129	gSystem->FreeDirectory(dir);
	130
586f2bc3	131	fAOD = 0;
	132	fChain->SetBranchAddress("AOD",&fAOD);
	133
	134	int nMax = fChain->GetEntries();
	135
be6e5811	136	printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax);
586f2bc3	137
	138	// set number of events in header
	139	if (fReaderHeader->GetLastEvent() == -1)
	140	fReaderHeader->SetLastEvent(nMax);
	141	else {
	142	Int_t nUsr = fReaderHeader->GetLastEvent();
	143	fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
	144	}
	145	}
	146
	147	//____________________________________________________________________________
	148
	149	void AliJetAODReader::ConnectTree(TTree* tree, TObject* /data/) {
	150	// Connect the tree
	151	// For AOD reader it's needed only to set the number of events
	152	fChain = (TChain*) tree;
	153
	154	Int_t nMax = fChain->GetEntries();
be6e5811	155	printf("\n AliJetAODReader: Total number of events in chain= %5d \n", nMax);
586f2bc3	156	// set number of events in header
	157	if (fReaderHeader->GetLastEvent() == -1)
	158	fReaderHeader->SetLastEvent(nMax);
	159	else {
	160	Int_t nUsr = fReaderHeader->GetLastEvent();
	161	fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
	162	}
	163	}
	164
	165	//____________________________________________________________________________
	166
9e4cc50d	167	Bool_t AliJetAODReader::FillMomentumArray()
586f2bc3	168	{
	169	// Clear momentum array
	170	ClearArray();
e53baffe	171	fRef->Clear();
921b3f4d	172
586f2bc3	173	fDebug = fReaderHeader->GetDebug();
	174
	175	if (!fAOD) {
	176	return kFALSE;
	177	}
	178
	179	// get number of tracks in event (for the loop)
	180	Int_t nt = fAOD->GetNTracks();
be6e5811	181	printf("AOD tracks: %5d \t", nt);
586f2bc3	182
	183	// temporary storage of signal and pt cut flag
	184	Int_t* sflag = new Int_t[nt];
	185	Int_t* cflag = new Int_t[nt];
	186
	187	// get cuts set by user
	188	Float_t ptMin = fReaderHeader->GetPtCut();
	189	Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
	190	Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
d61f057d	191	UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
e53baffe	192
586f2bc3	193	//loop over tracks
	194	Int_t aodTrack = 0;
	195	Float_t pt, eta;
	196	TVector3 p3;
d61f057d	197
586f2bc3	198	for (Int_t it = 0; it < nt; it++) {
d61f057d	199	AliAODTrack *track = fAOD->GetTrack(it);
e53baffe	200	UInt_t status = track->GetStatus();
d61f057d	201
	202	Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
	203	p3.SetXYZ(mom[0],mom[1],mom[2]);
	204	pt = p3.Pt();
	205	eta = p3.Eta();
e53baffe	206	if (status == 0) continue;
d61f057d	207	if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
d61f057d	208	if ( (eta > etaMax) \|\| (eta < etaMin)) continue; // checking eta cut
e53baffe	209
921b3f4d	210	if (aodTrack == 0) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
921b3f4d	211
e53baffe	212	new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
	213	sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
	214	cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
	215	aodTrack++;
d61f057d	216	fRef->Add(track);
586f2bc3	217	}
be6e5811	218	printf("Used AOD tracks: %5d \n", aodTrack);
586f2bc3	219	// set the signal flags
	220	fSignalFlag.Set(aodTrack,sflag);
	221	fCutFlag.Set(aodTrack,cflag);
d61f057d	222
680ed75f	223	delete [] sflag;
680ed75f	224	delete [] cflag;
d61f057d	225
586f2bc3	226	return kTRUE;
586f2bc3	227	}
be6e5811	228
	229	//__________________________________________________________
	230	void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
	231	{
	232	//
	233	// Set flag to apply MIP correction fApplyMIPCorrection
	234	// - exclusive with fApplyFractionHadronicCorrection
	235	//
	236
	237	fApplyMIPCorrection = val;
	238	if(fApplyMIPCorrection == kTRUE)
	239	{
	240	SetApplyFractionHadronicCorrection(kFALSE);
	241	printf("Enabling MIP Correction \n");
	242	}
	243	else
	244	{
	245	printf("Disabling MIP Correction \n");
	246	}
	247	}
	248
	249	//__________________________________________________________
	250	void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
	251	{
	252	//
	253	// Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
	254	// - exclusive with fApplyMIPCorrection
	255	//
	256
	257	fApplyFractionHadronicCorrection = val;
	258	if(fApplyFractionHadronicCorrection == kTRUE)
	259	{
	260	SetApplyMIPCorrection(kFALSE);
	261	printf("Enabling Fraction Hadronic Correction \n");
	262	}
	263	else
	264	{
	265	printf("Disabling Fraction Hadronic Correction \n");
	266	}
	267	}
	268
	269	//__________________________________________________________
	270	void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
	271	{
	272	//
	273	// Set value to fFractionHadronicCorrection (default is 0.3)
	274	// apply EMC hadronic correction fApplyFractionHadronicCorrection
	275	// - exclusive with fApplyMIPCorrection
	276	//
	277
	278	fFractionHadronicCorrection = val;
	279	if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
	280	{
	281	SetApplyFractionHadronicCorrection(kTRUE);
	282	printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
	283	}
	284	else
	285	{
	286	SetApplyFractionHadronicCorrection(kFALSE);
	287	}
	288	}
	289
	290	//____________________________________________________________________________
	291	void AliJetAODReader::CreateTasks(TChain* tree)
292	{
4751efb5	293	//
	294	// For reader task initialization
	295	//
	296
be6e5811	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();
f23f7908	452	Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
be6e5811	453	if(fArrayInitialised) fUnitArray->Delete();
	454
	455	if(fTpcGrid->GetGridType()==0)
	456	{ // Fill the following quantities :
f23f7908	457	// Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
be6e5811	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);
f23f7908	464	deltaEta = fTpcGrid->GetDeta();
	465	deltaPhi = fTpcGrid->GetDphi();
	466	new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	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);
f23f7908	522	deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
	523	deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
	524	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	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);
f23f7908	530	deltaEta = fTpcGrid->GetDeta();
	531	deltaPhi = fTpcGrid->GetDphi();
	532	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	533	}
	534	else {
	535	if(fDZ) {
	536	if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
	537	if(nBin<fNumUnits+nElements+n0)
	538	{
f23f7908	539	phi = eta = 0.;
be6e5811	540	fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
f23f7908	541	deltaEta = fGrid0->GetDeta();
	542	deltaPhi = fGrid0->GetDphi();
	543	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	544	}
	545	else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
	546	{
f23f7908	547	phi = eta = 0.;
be6e5811	548	fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
f23f7908	549	deltaEta = fGrid1->GetDeta();
	550	deltaPhi = fGrid1->GetDphi();
	551	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	552	}
	553	else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
	554	{
f23f7908	555	phi = eta = 0.;
be6e5811	556	fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
f23f7908	557	deltaEta = fGrid2->GetDeta();
	558	deltaPhi = fGrid2->GetDphi();
	559	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	560	}
	561	else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
	562	{
f23f7908	563	phi = eta = 0.;
be6e5811	564	fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
f23f7908	565	deltaEta = fGrid3->GetDeta();
	566	deltaPhi = fGrid3->GetDphi();
	567	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	568	}
	569	else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
	570	{
f23f7908	571	phi = eta = 0.;
be6e5811	572	fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
f23f7908	573	deltaEta = fGrid4->GetDeta();
	574	deltaPhi = fGrid4->GetDphi();
	575	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	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	}