[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

    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., 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();
586f2bc3	172	fDebug = fReaderHeader->GetDebug();
	173
	174	if (!fAOD) {
	175	return kFALSE;
	176	}
	177
	178	// get number of tracks in event (for the loop)
	179	Int_t nt = fAOD->GetNTracks();
be6e5811	180	printf("AOD tracks: %5d \t", nt);
586f2bc3	181
	182	// temporary storage of signal and pt cut flag
	183	Int_t* sflag = new Int_t[nt];
	184	Int_t* cflag = new Int_t[nt];
	185
	186	// get cuts set by user
	187	Float_t ptMin = fReaderHeader->GetPtCut();
	188	Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
	189	Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
d61f057d	190	UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
e53baffe	191
586f2bc3	192	//loop over tracks
	193	Int_t aodTrack = 0;
	194	Float_t pt, eta;
	195	TVector3 p3;
d61f057d	196
586f2bc3	197	for (Int_t it = 0; it < nt; it++) {
d61f057d	198	AliAODTrack *track = fAOD->GetTrack(it);
e53baffe	199	UInt_t status = track->GetStatus();
d61f057d	200
	201	Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
	202	p3.SetXYZ(mom[0],mom[1],mom[2]);
	203	pt = p3.Pt();
	204	eta = p3.Eta();
e53baffe	205	if (status == 0) continue;
d61f057d	206	if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
d61f057d	207	if ( (eta > etaMax) \|\| (eta < etaMin)) continue; // checking eta cut
e53baffe	208
	209	new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
	210	sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
	211	cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
	212	aodTrack++;
d61f057d	213	fRef->Add(track);
586f2bc3	214	}
be6e5811	215	printf("Used AOD tracks: %5d \n", aodTrack);
586f2bc3	216	// set the signal flags
	217	fSignalFlag.Set(aodTrack,sflag);
	218	fCutFlag.Set(aodTrack,cflag);
d61f057d	219
680ed75f	220	delete [] sflag;
680ed75f	221	delete [] cflag;
d61f057d	222
586f2bc3	223	return kTRUE;
586f2bc3	224	}
be6e5811	225
	226	//__________________________________________________________
	227	void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
	228	{
	229	//
	230	// Set flag to apply MIP correction fApplyMIPCorrection
	231	// - exclusive with fApplyFractionHadronicCorrection
	232	//
	233
	234	fApplyMIPCorrection = val;
	235	if(fApplyMIPCorrection == kTRUE)
	236	{
	237	SetApplyFractionHadronicCorrection(kFALSE);
	238	printf("Enabling MIP Correction \n");
	239	}
	240	else
	241	{
	242	printf("Disabling MIP Correction \n");
	243	}
	244	}
	245
	246	//__________________________________________________________
	247	void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
	248	{
	249	//
	250	// Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
	251	// - exclusive with fApplyMIPCorrection
	252	//
	253
	254	fApplyFractionHadronicCorrection = val;
	255	if(fApplyFractionHadronicCorrection == kTRUE)
	256	{
	257	SetApplyMIPCorrection(kFALSE);
	258	printf("Enabling Fraction Hadronic Correction \n");
	259	}
	260	else
	261	{
	262	printf("Disabling Fraction Hadronic Correction \n");
	263	}
	264	}
	265
	266	//__________________________________________________________
	267	void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
	268	{
	269	//
	270	// Set value to fFractionHadronicCorrection (default is 0.3)
	271	// apply EMC hadronic correction fApplyFractionHadronicCorrection
	272	// - exclusive with fApplyMIPCorrection
	273	//
	274
	275	fFractionHadronicCorrection = val;
	276	if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
	277	{
	278	SetApplyFractionHadronicCorrection(kTRUE);
	279	printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
	280	}
	281	else
	282	{
	283	SetApplyFractionHadronicCorrection(kFALSE);
	284	}
	285	}
	286
	287	//____________________________________________________________________________
	288	void AliJetAODReader::CreateTasks(TChain* tree)
289	{
290	fDebug = fReaderHeader->GetDebug();
291	fDZ = fReaderHeader->GetDZ();
292	fTree = tree;
293
294	// Init EMCAL geometry and create UnitArray object
295	SetEMCALGeometry();
296	// cout << "In create task" << endl;
297	InitParameters();
298	InitUnitArray();
299
300	fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
301	fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
302	fFillUAFromTracks->SetReaderHeader(fReaderHeader);
303	fFillUAFromTracks->SetGeom(fGeom);
304	fFillUAFromTracks->SetTPCGrid(fTpcGrid);
305	fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
306
307	if(fDZ)
308	{
309	fFillUAFromTracks->SetGrid0(fGrid0);
310	fFillUAFromTracks->SetGrid1(fGrid1);
311	fFillUAFromTracks->SetGrid2(fGrid2);
312	fFillUAFromTracks->SetGrid3(fGrid3);
313	fFillUAFromTracks->SetGrid4(fGrid4);
314	}
315	fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
316	fFillUAFromTracks->SetHadCorrector(fHadCorr);
317	fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
318	fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
319	fFillUAFromEMCalDigits->SetGeom(fGeom);
320	fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
321	fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
322	fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
323	fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
324	fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
325
326	fFillUnitArray->Add(fFillUAFromTracks);
327	fFillUnitArray->Add(fFillUAFromEMCalDigits);
328	fFillUAFromTracks->SetActive(kFALSE);
329	fFillUAFromEMCalDigits->SetActive(kFALSE);
330
331	cout << "Tasks instantiated at that stage ! " << endl;
332	cout << "You can loop over events now ! " << endl;
333
334	}
335
336	//____________________________________________________________________________
337	Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
338	{
339	//
340	// Main function
341	// Fill the reader part
342	//
343
344	fProcId = procid;
345	fRefArray = refArray;
346	//(not used ?) Int_t nEntRef = fRefArray->GetEntries();
347	//(not used ?) Int_t nEntUnit = fUnitArray->GetEntries();
348
349	// clear momentum array
350	ClearArray();
351
352	fDebug = fReaderHeader->GetDebug();
353	fOpt = fReaderHeader->GetDetector();
354
355	if(!fAOD) {
356	return kFALSE;
357	}
358
359	// TPC only or Digits+TPC or Clusters+TPC
360	if(fOpt%2==!0 && fOpt!=0){
361	fFillUAFromTracks->SetAOD(fAOD);
362	fFillUAFromTracks->SetActive(kTRUE);
363	fFillUAFromTracks->SetUnitArray(fUnitArray);
364	fFillUAFromTracks->SetRefArray(fRefArray);
365	fFillUAFromTracks->SetProcId(fProcId);
366	// fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
367	fFillUAFromTracks->Exec("tpc");
368	if(fOpt==1){
369	fNumCandidate = fFillUAFromTracks->GetMult();
370	fNumCandidateCut = fFillUAFromTracks->GetMultCut();
371	}
372	}
373
374	// Digits only or Digits+TPC
375	if(fOpt>=2 && fOpt<=3){
376	fFillUAFromEMCalDigits->SetAOD(fAOD);
377	fFillUAFromEMCalDigits->SetActive(kTRUE);
378	fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
379	fFillUAFromEMCalDigits->SetRefArray(fRefArray);
380	fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
381	fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
382	fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
383	fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
384	fNumCandidate = fFillUAFromEMCalDigits->GetMult();
385	fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
386	}
387
388	// fFillUnitArray->ExecuteTask(); // => Temporarily commented
389
390	return kTRUE;
391	}
392
393	//____________________________________________________________________________
394	Bool_t AliJetAODReader::SetEMCALGeometry()
395	{
396	//
397	// Set the EMCal Geometry
398	//
399
400	if (!fTree->GetFile())
401	return kFALSE;
402
403	TString geomFile(fTree->GetFile()->GetName());
404	geomFile.ReplaceAll("AliESDs", "geometry");
405
406	// temporary workaround for PROOF bug #18505
407	geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
408	if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
409
410	// Define EMCAL geometry to be able to read ESDs
411	fGeom = AliJetDummyGeo::GetInstance();
412	if (fGeom == 0)
413	fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
414
415	// To be setted to run some AliEMCALGeometry functions
416	TGeoManager::Import(geomFile);
417	fGeom->GetTransformationForSM();
418	printf("\n EMCal Geometry set ! \n");
419
420	return kTRUE;
421
422	}
423
424	//____________________________________________________________________________
425	void AliJetAODReader::InitParameters()
426	{
427	// Initialise parameters
428	fOpt = fReaderHeader->GetDetector();
429	// fHCorrection = 0; // For hadron correction
430	fHadCorr = 0; // For hadron correction
431	if(fEFlag==kFALSE){
432	if(fOpt==0 \|\| fOpt==1)
433	fECorrection = 0; // For electron correction
434	else fECorrection = 1; // For electron correction
435	}
436	fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL
437	if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
438	}
439
440	//____________________________________________________________________________
441	void AliJetAODReader::InitUnitArray()
442	{
443	//Initialises unit arrays
444	Int_t nElements = fTpcGrid->GetNEntries();
f23f7908	445	Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
be6e5811	446	if(fArrayInitialised) fUnitArray->Delete();
	447
	448	if(fTpcGrid->GetGridType()==0)
	449	{ // Fill the following quantities :
f23f7908	450	// Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
be6e5811	451	// detector flag, in/out jet, pt cut, mass, cluster ID)
	452	for(Int_t nBin = 1; nBin < nElements+1; nBin++)
	453	{
	454	// fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
	455	fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
	456	phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
f23f7908	457	deltaEta = fTpcGrid->GetDeta();
	458	deltaPhi = fTpcGrid->GetDphi();
	459	new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	460	}
	461	}
	462
	463	if(fTpcGrid->GetGridType()==1)
	464	{
	465	Int_t nGaps = 0;
	466	Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
	467
	468	if(fDZ)
	469	{
	470	// Define a grid of cell for the gaps between SM
	471	Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
	472	Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
	473	fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
	474	fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
	475	fGrid0->SetGridType(0);
	476	fGrid0->SetMatrixIndexes();
	477	fGrid0->SetIndexIJ();
	478	n0 = fGrid0->GetNEntries();
	479	fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
	480	fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
	481	fGrid1->SetGridType(0);
	482	fGrid1->SetMatrixIndexes();
	483	fGrid1->SetIndexIJ();
	484	n1 = fGrid1->GetNEntries();
	485	fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
	486	fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
	487	fGrid2->SetGridType(0);
	488	fGrid2->SetMatrixIndexes();
	489	fGrid2->SetIndexIJ();
	490	n2 = fGrid2->GetNEntries();
	491	fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
	492	fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
	493	fGrid3->SetGridType(0);
	494	fGrid3->SetMatrixIndexes();
	495	fGrid3->SetIndexIJ();
	496	n3 = fGrid3->GetNEntries();
	497	fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
	498	fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
	499	fGrid4->SetGridType(0);
	500	fGrid4->SetMatrixIndexes();
	501	fGrid4->SetIndexIJ();
	502	n4 = fGrid4->GetNEntries();
	503
	504	nGaps = n0+n1+n2+n3+n4;
	505
	506	}
	507
	508	for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
	509	{
	510	if(nBin<fNumUnits)
	511	{
	512	fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
	513	// fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
	514	phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
f23f7908	515	deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
	516	deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
	517	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	518	}
	519	else {
	520	if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
	521	fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
	522	phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
f23f7908	523	deltaEta = fTpcGrid->GetDeta();
	524	deltaPhi = fTpcGrid->GetDphi();
	525	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	526	}
	527	else {
	528	if(fDZ) {
	529	if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
	530	if(nBin<fNumUnits+nElements+n0)
	531	{
f23f7908	532	phi = eta = 0.;
be6e5811	533	fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
f23f7908	534	deltaEta = fGrid0->GetDeta();
	535	deltaPhi = fGrid0->GetDphi();
	536	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	537	}
	538	else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
	539	{
f23f7908	540	phi = eta = 0.;
be6e5811	541	fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
f23f7908	542	deltaEta = fGrid1->GetDeta();
	543	deltaPhi = fGrid1->GetDphi();
	544	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	545	}
	546	else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
	547	{
f23f7908	548	phi = eta = 0.;
be6e5811	549	fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
f23f7908	550	deltaEta = fGrid2->GetDeta();
	551	deltaPhi = fGrid2->GetDphi();
	552	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	553	}
	554	else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
	555	{
f23f7908	556	phi = eta = 0.;
be6e5811	557	fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
f23f7908	558	deltaEta = fGrid3->GetDeta();
	559	deltaPhi = fGrid3->GetDphi();
	560	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	561	}
	562	else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
	563	{
f23f7908	564	phi = eta = 0.;
be6e5811	565	fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
f23f7908	566	deltaEta = fGrid4->GetDeta();
	567	deltaPhi = fGrid4->GetDphi();
	568	new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811	569	}
	570	}
	571	} // end if(fDZ)
	572	} // end else 2
	573	} // end else 1
	574	} // end loop on nBin
	575	} // end grid type == 1
	576	fArrayInitialised = 1;
	577	}