[u/mrichter/AliRoot.git] / JETAN / AliJetESDReader.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//------------------------------------------------------------------------- 
// Jet ESD Reader 
// ESD reader for jet analysis
// Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
//          Magali Estienne <magali.estienne@subatech.in2p3.fr>
//------------------------------------------------------------------------- 

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

// --- ROOT system ---
#include <TSystem.h>
#include <TStopwatch.h>
#include <TLorentzVector.h>
#include <TVector3.h>
#include "TTask.h"
#include "TTree.h"
#include "TFile.h"
#include <TGeoManager.h>
#include <assert.h>
#include <TRefArray.h>
#include <TMath.h>
#include <TProcessID.h>
#include <TRandom3.h>

// --- AliRoot header files ---
#include "AliJetESDReader.h"
#include "AliJetESDReaderHeader.h"
#include "AliESDEvent.h"
#include "AliESD.h"
#include "AliESDtrack.h"
#include "AliJetDummyGeo.h"
#include "AliJetFillUnitArrayTracks.h"
#include "AliJetFillUnitArrayEMCalDigits.h"
#include "AliJetUnitArray.h"
#include "AliAnalysisTask.h"

ClassImp(AliJetESDReader)

AliJetESDReader::AliJetESDReader():
  AliJetReader(),  
  fGeom(0),
  fChain(0x0),
  fTree(0x0),
  fESD(0x0),
  fHadCorr(0x0),
  fTpcGrid(0x0),
  fEmcalGrid(0x0),
  fGrid0(0),
  fGrid1(0),
  fGrid2(0),
  fGrid3(0),
  fGrid4(0),
  fPtCut(0),
  fHCorrection(0),
  fECorrection(0),
  fEFlag(kFALSE),
  fNumUnits(0),
  fDebug(0),
  fMass(0),
  fSign(0),
  fNIn(0),
  fOpt(0),
  fDZ(0),
  fNeta(0),
  fNphi(0),
  fArrayInitialised(0),
  fRefArray(0x0),
  fProcId(kFALSE)
{
  // Constructor 
}

//____________________________________________________________________________
AliJetESDReader::~AliJetESDReader()
{
  // Destructor
    delete fChain;
    delete fTree;
    delete fESD;
    delete fTpcGrid;
    delete fEmcalGrid;
    if(fDZ)
      {
	delete fGrid0;
	delete fGrid1;
	delete fGrid2;
	delete fGrid3;
	delete fGrid4;
      }
}

//____________________________________________________________________________
void AliJetESDReader::OpenInputFiles()
{
  // chain for the ESDs
  fChain   = new TChain("esdTree");

  // get directory and pattern name from the header
   const char* dirName=fReaderHeader->GetDirectory();
   const char* pattern=fReaderHeader->GetPattern();
    
   // Add files matching patters to the chain
   void *dir  = gSystem->OpenDirectory(dirName);
   const char *name = 0x0;
   int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
   int a = 0;
   while ((name = gSystem->GetDirEntry(dir))){
       if (a>=nesd) continue;
       
       if (strstr(name,pattern)){
	 printf("Adding %s\n",name);
	 char path[256];
	 sprintf(path,"%s/%s/AliESDs.root",dirName,name);
	 fChain->AddFile(path);
	 a++;
       }
   }
  
  gSystem->FreeDirectory(dir);
  
  fESD = new AliESDEvent();
  fESD->ReadFromTree(fChain);

  int nMax = fChain->GetEntries(); 

  printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);

  // set number of events in header
  if (fReaderHeader->GetLastEvent() == -1)
    fReaderHeader->SetLastEvent(nMax);
  else {
    Int_t nUsr = fReaderHeader->GetLastEvent();
    fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
  }

}

//____________________________________________________________________________
void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) {
    // Connect the tree
     fESD   = (AliESDEvent*) esd;
}

//____________________________________________________________________________
Bool_t AliJetESDReader::FillMomentumArray()
{
  // Fill momentum array

  Int_t goodTrack = 0;
  Int_t nt = 0;
  Float_t pt, eta;
  TVector3 p3;
  
  // clear momentum array
  ClearArray();
  fDebug = fReaderHeader->GetDebug();
  fOpt = fReaderHeader->GetDetector();

  if (!fESD) {
      return kFALSE;
  }

  // get number of tracks in event (for the loop)
  nt = fESD->GetNumberOfTracks();
  printf("Fill Momentum Array %5d  \n", nt);
 
  // temporary storage of signal and pt cut flag
  Int_t* sflag  = new Int_t[nt];
  Int_t* cflag  = new Int_t[nt];
  
  // get cuts set by user
  Float_t ptMin =  fReaderHeader->GetPtCut();
  Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
  Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
  
  //loop over tracks in ESD
  for (Int_t it = 0; it < nt; it++) {
      AliESDtrack *track = fESD->GetTrack(it);
      UInt_t status = track->GetStatus();
      
      Double_t mom[3];
      track->GetPxPyPz(mom);
      p3.SetXYZ(mom[0],mom[1],mom[2]);
      pt = p3.Pt();
      if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
      if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
      if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() 
	  && TMath::Abs(track->GetLabel()) > 10000)  continue;      // quality check
      if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() 
	  && TMath::Abs(track->GetLabel()) < 10000)  continue;      // quality check
      eta = p3.Eta();
      if ( (eta > etaMax) || (eta < etaMin))         continue;      // checking eta cut
      
      new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
      sflag[goodTrack]=0;
      if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
      cflag[goodTrack]=0;
      if (pt > ptMin) cflag[goodTrack]=1;                           // pt cut
      goodTrack++;
  }

  // set the signal flags
  fSignalFlag.Set(goodTrack,sflag);
  fCutFlag.Set(goodTrack,cflag);

  delete[] sflag;
  delete[] cflag;

  return kTRUE;

}

//____________________________________________________________________________
void AliJetESDReader::CreateTasks(TChain* tree)
{
  //
  // For reader task initialization
  //

  fDebug = fReaderHeader->GetDebug();
  fDZ = fReaderHeader->GetDZ();
  fTree = tree;

  // Init EMCAL geometry 
  SetEMCALGeometry();
  // Init parameters
  InitParameters();
  // Create and init unit array
  InitUnitArray();

  // Create global reader task for analysis 
  fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
  // Create a task for to fill the charged particle information 
  fFillUAFromTracks = new AliJetFillUnitArrayTracks(); 
  fFillUAFromTracks->SetReaderHeader(fReaderHeader);
  fFillUAFromTracks->SetGeom(fGeom);
  fFillUAFromTracks->SetTPCGrid(fTpcGrid);
  fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
  if(fDZ)
    { // Calo dead zones inclusion
      fFillUAFromTracks->SetGrid0(fGrid0);
      fFillUAFromTracks->SetGrid1(fGrid1);
      fFillUAFromTracks->SetGrid2(fGrid2);
      fFillUAFromTracks->SetGrid3(fGrid3);
      fFillUAFromTracks->SetGrid4(fGrid4);
    }
  fFillUAFromTracks->SetHadCorrection(fHCorrection);
  fFillUAFromTracks->SetHadCorrector(fHadCorr);
  // Create a task for to fill the neutral particle information 
  fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits();
  fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
  fFillUAFromEMCalDigits->SetGeom(fGeom);
  fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
  fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
  fFillUAFromEMCalDigits->SetEleCorrection(fECorrection);
  // Add the task to global task
  fFillUnitArray->Add(fFillUAFromTracks);
  fFillUnitArray->Add(fFillUAFromEMCalDigits);
  fFillUAFromTracks->SetActive(kFALSE);
  fFillUAFromEMCalDigits->SetActive(kFALSE);

  cout << "Tasks instantiated at that stage ! " << endl;
  cout << "You can loop over events now ! " << endl;
   
}

//____________________________________________________________________________
Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
{
  //
  // Reader task execussion
  //

  fProcId = procid;
  fRefArray = refArray;
  vector<Float_t> vtmp(3);

  // clear momentum array
  ClearArray();

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

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

  // Digits only or Digits+TPC  
  if(fOpt>=2 && fOpt<=3){
    fFillUAFromEMCalDigits->SetESD(fESD);
    fFillUAFromEMCalDigits->SetActive(kTRUE);
    fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
    fFillUAFromEMCalDigits->SetRefArray(fRefArray);
    fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
    fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
    fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
    fFillUAFromEMCalDigits->Exec("digits");      // Temporarily added
    fNumCandidate = fFillUAFromEMCalDigits->GetMult();
    fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
  }

  //  fFillUnitArray->ExecuteTask();             // Temporarily commented

  return kTRUE;
}

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

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

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

  // Define EMCAL geometry to be able to read ESDs
  fGeom = AliJetDummyGeo::GetInstance();
  if (fGeom == 0)
    fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
  
  // To be setted to run some AliEMCALGeometry functions
  TGeoManager::Import(geomFile);
  fGeom->GetTransformationForSM();  
  printf("\n EMCal Geometry set ! \n");
  
  return kTRUE;
}

//____________________________________________________________________________  
void AliJetESDReader::InitParameters()
{
  // Initialise parameters
  fOpt = fReaderHeader->GetDetector();
  fHadCorr        = 0;          // For hadron correction
  if(fEFlag==kFALSE){
    if(fOpt==0 || fOpt==1)  
      fECorrection    = 0;        // For electron correction
    else fECorrection = 1;        // For electron correction
  }
  fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
  if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
}

//____________________________________________________________________________
void AliJetESDReader::InitUnitArray()
{
  //
  // Create and Initialise unit arrays
  //

  Int_t nElements = fTpcGrid->GetNEntries();
  Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
  if(fArrayInitialised) fUnitArray->Delete();

  if(fTpcGrid->GetGridType()==0)
    { // Fill the following quantities :
      // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi, 
      // detector flag, in/out jet, pt cut, mass, cluster ID)
      for(Int_t nBin = 1; nBin < nElements+1; nBin++)
	{
	  fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
	  phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
	  Deta = fTpcGrid->GetDeta();
	  Dphi = fTpcGrid->GetDphi();
	  new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
	}
    }

  if(fTpcGrid->GetGridType()==1)
    {
      Int_t nGaps = 0;
      Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;

      if(fDZ)
	{
	  // Define a grid of cell for the gaps between SM
	  Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
	  Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
	  fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
	  fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
	  fGrid0->SetGridType(0);
	  fGrid0->SetMatrixIndexes();
	  fGrid0->SetIndexIJ();
	  n0 = fGrid0->GetNEntries();
	  fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
	  fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
	  fGrid1->SetGridType(0);
	  fGrid1->SetMatrixIndexes();
	  fGrid1->SetIndexIJ();
	  n1 = fGrid1->GetNEntries();
	  fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
	  fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
	  fGrid2->SetGridType(0);
	  fGrid2->SetMatrixIndexes();
	  fGrid2->SetIndexIJ();
	  n2 = fGrid2->GetNEntries();
	  fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
	  fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
	  fGrid3->SetGridType(0);  
	  fGrid3->SetMatrixIndexes();
	  fGrid3->SetIndexIJ();
	  n3 = fGrid3->GetNEntries();
	  fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
	  fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
	  fGrid4->SetGridType(0);
	  fGrid4->SetMatrixIndexes();
	  fGrid4->SetIndexIJ();
	  n4 = fGrid4->GetNEntries();

	  if(fDebug>1) 
	    {
	      cout << "n0 cells: " << n0 << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl;
	      cout << "n1 cells: " << n1 << "phimin1: " << phimin1 << ", phimax1: " << phimax1 << endl;
	      cout << "n2 cells: " << n2 << "phimin2: " << phimin2 << ", phimax2: " << phimax2 << endl;
	      cout << "n3 cells: " << n3 << "phimin3: " << phimin3 << ", phimax3: " << phimax3 << endl;
	      cout << "n4 cells: " << n4 << "phimin4: " << phimin4 << ", phimax4: " << phimax4 << endl;
	    }
	  
	  nGaps = n0+n1+n2+n3+n4;

	}

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