Updates provided by Magali Estienne
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Jul 2009 09:10:18 +0000 (09:10 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Jul 2009 09:10:18 +0000 (09:10 +0000)
1 - remove the coding convention violations you sent to me few weeks ago
2 - look for jets using charged + neutral information from AODs

26 files changed:
JETAN/AliJet.cxx
JETAN/AliJet.h
JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx [new file with mode: 0755]
JETAN/AliJetAODFillUnitArrayEMCalDigits.h [new file with mode: 0755]
JETAN/AliJetAODFillUnitArrayTracks.cxx [new file with mode: 0644]
JETAN/AliJetAODFillUnitArrayTracks.h [new file with mode: 0644]
JETAN/AliJetAODReader.cxx
JETAN/AliJetAODReader.h
JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx [new file with mode: 0755]
JETAN/AliJetESDFillUnitArrayEMCalDigits.h [new file with mode: 0755]
JETAN/AliJetESDFillUnitArrayTracks.cxx [new file with mode: 0644]
JETAN/AliJetESDFillUnitArrayTracks.h [new file with mode: 0644]
JETAN/AliJetESDReader.cxx
JETAN/AliJetESDReader.h
JETAN/AliJetFillUnitArray.cxx [new file with mode: 0644]
JETAN/AliJetFillUnitArray.h [new file with mode: 0644]
JETAN/AliJetFinder.cxx
JETAN/AliJetFinder.h
JETAN/AliJetReader.cxx
JETAN/AliJetReader.h
JETAN/AliJetUnitArray.cxx
JETAN/AliJetUnitArray.h
JETAN/AliUA1JetFinderV2.cxx
JETAN/AliUA1JetFinderV2.h
JETAN/JETANLinkDef.h
JETAN/libJETAN.pkg

index 23e4e68..4367b66 100644 (file)
@@ -23,7 +23,6 @@
 
  
 #include <Riostream.h>
-//#include <vector.h>
 #include <TClonesArray.h>
 #include <TLorentzVector.h>
 
@@ -44,12 +43,8 @@ AliJet::AliJet():
   fPtIn(0),
   fPtChPtCutIn(0),
   fEnTotChPtCutIn(0),
-  fVectorSizeIn(0),
   fDetIn(0),
-  fVPx(0),
-  fVPy(0),
-  fVPz(0)
-             //  fVectorIn(0)
+  fTrackRef(new TRefArray())
 {
   // Default constructor
   fJets = new TClonesArray("TLorentzVector",1000);
@@ -62,7 +57,6 @@ AliJet::AliJet():
   fNCells = TArrayI();
   fPtChPtCutIn = TArrayF();
   fEnTotChPtCutIn = TArrayF();
-  fVectorSizeIn = TArrayI();
   fDetIn = TArrayI();
 } 
 
@@ -280,34 +274,6 @@ void AliJet::SetDetectorFlagIn(Int_t* x)
 
 ////////////////////////////////////////////////////////////////////////
 
-void AliJet::SetVectorSizeIn(Int_t* x)
-{
-  if (fNInput>0) fVectorSizeIn.Set(fNInput, x);
-}
-
-////////////////////////////////////////////////////////////////////////
-
-void AliJet::SetVectorPxIn(vector< vector<Float_t> > pxT)
-{
-  fVPx = pxT;;
-}
-
-////////////////////////////////////////////////////////////////////////
-
-void AliJet::SetVectorPyIn(vector< vector<Float_t> > pyT)
-{
-  fVPy = pyT;;
-}
-
-////////////////////////////////////////////////////////////////////////
-
-void AliJet::SetVectorPzIn(vector< vector<Float_t> > pzT)
-{
-  fVPz = pzT;;
-}
-
-////////////////////////////////////////////////////////////////////////
-
 void AliJet::SetPtFromSignal(Float_t* p)
 {
   // set information of percentage of pt of jets
@@ -348,7 +314,6 @@ void AliJet::ClearJets(Option_t *option)
   fNCells.Set(0);
   fPtChPtCutIn.Set(0);
   fEnTotChPtCutIn.Set(0);
-  fVectorSizeIn.Set(0);
   fDetIn.Set(0);
 }
 
index 74689e9..2d05825 100644 (file)
 //---------------------------------------------------------------------
 
 #include <Riostream.h>
-#include <vector>
  
 #include <TObject.h>
 #include <TArrayI.h>
 #include <TArrayF.h>
+#include <TRefArray.h>
 
 class TClonesArray;
 class TLorentzVector;
@@ -41,10 +41,7 @@ class AliJet : public TObject
   TArrayF       GetPtIn() const { return fPtIn; }
   TArrayF       GetPtChargedPtCutIn() const { return fPtChPtCutIn; }
   TArrayF       GetEnTotChargedPtCutIn() const {return fEnTotChPtCutIn; }
-  TArrayI       GetVectorSizeIn() const { return fVectorSizeIn;  }
-  vector< vector<Float_t> > GetVectorPxIn() const { return fVPx; }
-  vector< vector<Float_t> > GetVectorPyIn() const { return fVPy; }
-  vector< vector<Float_t> > GetVectorPzIn() const { return fVPz; }
+  TRefArray*    GetTrackRef() const { return fTrackRef;}
   TArrayI       GetDetectorFlagIn() const { return fDetIn; } 
   Double_t      GetEtAvg() const { return fEtAvg; }
 
@@ -74,10 +71,7 @@ class AliJet : public TObject
   void SetInJet(Int_t* idx);
   void SetPtChargedPtCutIn(Float_t* pt2T);
   void SetEnTotChargedPtCutIn(Float_t* en2T);
-  void SetVectorSizeIn(Int_t* vectT);
-  void SetVectorPxIn(vector< vector<Float_t> > pxT);
-  void SetVectorPyIn(vector< vector<Float_t> > pyT);
-  void SetVectorPzIn(vector< vector<Float_t> > pzT);
+  void SetTrackReferences(TRefArray* ref) {fTrackRef = ref;}
   void SetDetectorFlagIn(Int_t* detT);
   void SetEtAvg(Double_t et) { fEtAvg = et; }
   // others
@@ -105,11 +99,8 @@ class AliJet : public TObject
   TArrayF  fPtIn;                 // Arrays of input particles kine:Pt
   TArrayF  fPtChPtCutIn;          // Arrays of input particles kin:Pt Charged with pt cut
   TArrayF  fEnTotChPtCutIn;       // Arrays of total energy with pt cut on charged + cut min on cell
-  TArrayI  fVectorSizeIn;         // Arrays of number of charged tracks in each unitArray
   TArrayI  fDetIn;                // Arrays of detector type of each UnitArray
-  vector< vector<Float_t> > fVPx; //||
-  vector< vector<Float_t> > fVPy; //||
-  vector< vector<Float_t> > fVPz; //||
+  TRefArray* fTrackRef;           //|| Reference to tracks which could belongs to the jet
 
   ClassDef(AliJet,2)
 };
diff --git a/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx b/JETAN/AliJetAODFillUnitArrayEMCalDigits.cxx
new file mode 100755 (executable)
index 0000000..2e9429d
--- /dev/null
@@ -0,0 +1,390 @@
+
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+//======================================================================
+// Fill Unit Array 
+// Called by AOD reader for jet analysis
+// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
+//======================================================================
+
+class TSystem;
+class TLorentzVector;
+class TGeoManager;
+class TArrayS;
+
+// --- AliRoot header files ---
+#include "AliJetUnitArray.h"
+#include "AliJetAODFillUnitArrayEMCalDigits.h"
+class AliJetFinder;
+
+// #include "AliEMCALCalibData.h"
+// #include "AliCDBManager.h"
+// class AliCDBStorage;
+// #include "AliCDBEntry.h"
+
+
+ClassImp(AliJetAODFillUnitArrayEMCalDigits)
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits()
+  : AliJetFillUnitArray(),
+    fAOD(0),
+    fNIn(0),
+    fCluster(0),
+    fNCEMCAL(0),
+    fNCPHOS(0),
+    fNCCalo(0),
+    fApplyElectronCorrection(kFALSE),
+    fApplyFractionHadronicCorrection(kFALSE),
+    fFractionHadronicCorrection(0.3),
+    fClus(0x0),
+    fNDigitEmcal(0),
+    fNDigitEmcalCut(0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(AliAODEvent *aod)
+  : AliJetFillUnitArray(),
+    fAOD(aod),
+    fNIn(0),
+    fCluster(0),
+    fNCEMCAL(0),
+    fNCPHOS(0),
+    fNCCalo(0),
+    fApplyElectronCorrection(kFALSE),
+    fApplyFractionHadronicCorrection(kFALSE),
+    fFractionHadronicCorrection(0.3),
+    fClus(0x0),
+    fNDigitEmcal(0),
+    fNDigitEmcalCut(0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(const AliJetAODFillUnitArrayEMCalDigits &det)
+  : AliJetFillUnitArray(det),
+    fAOD(det.fAOD),
+    fNIn(det.fNIn),
+    fCluster(det.fCluster),
+    fNCEMCAL(det.fNCEMCAL),
+    fNCPHOS(det.fNCPHOS),
+    fNCCalo(det.fNCCalo),
+    fApplyElectronCorrection(det.fApplyElectronCorrection),
+    fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
+    fFractionHadronicCorrection(det.fFractionHadronicCorrection),
+    fClus(det.fClus),
+    fNDigitEmcal(det.fNDigitEmcal),
+    fNDigitEmcalCut(det.fNDigitEmcalCut)
+{
+  // Copy constructor
+}
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayEMCalDigits& AliJetAODFillUnitArrayEMCalDigits::operator=(const AliJetAODFillUnitArrayEMCalDigits& other)
+{
+  // Assignment
+
+  fAOD = other.fAOD;
+  fNIn = other.fNIn;
+  fCluster = other.fCluster;
+  fNCEMCAL = other.fNCEMCAL;
+  fNCPHOS = other.fNCPHOS;
+  fNCCalo = other.fNCCalo;
+  fApplyElectronCorrection = other.fApplyElectronCorrection;
+  fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection;
+  fFractionHadronicCorrection = other.fFractionHadronicCorrection;
+  fClus = other.fClus;
+  fNDigitEmcal = other.fNDigitEmcal;
+  fNDigitEmcalCut = other.fNDigitEmcalCut;
+
+  return (*this);
+
+}
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayEMCalDigits::~AliJetAODFillUnitArrayEMCalDigits()
+{
+  // destructor
+}
+
+
+
+//_____________________________________________________________________________
+void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
+//void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* option)
+{
+  //
+  // Main method.
+  // Fill the unit array with the neutral particle information from the EMCal cells in AOD
+
+  fDebug = fReaderHeader->GetDebug();
+  fOpt = fReaderHeader->GetDetector();
+  fCluster = fReaderHeader->GetCluster();
+
+  // Init parameters
+  InitParameters();
+
+//(not used ?)  Int_t   nDigitTot      = 0;
+  Int_t   goodDigit      = 0;
+//(not used ?)  Int_t   beg            = 0;
+//(not used ?)  Int_t   end            = 0; 
+  Int_t   index          = 0;
+//(not used ?)  Int_t   count          = 0;
+
+//(not used ?)  Int_t nRefEnt = fRefArray->GetEntries();
+
+  if(!fCluster) { // Keep all digit information
+    // Loop over all cell information
+    //------------------------------------------------------------------
+    AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
+    Int_t ncell = cells.GetNumberOfCells() ;
+//(not used ?)    Int_t type = cells.GetType();
+
+    for (Int_t icell=  0; icell <  ncell; icell++) {
+      Int_t   digitID   = cells.GetCellNumber(icell);
+      Float_t digitAmp  = cells.GetAmplitude(icell);
+      Float_t digitEn   = digitAmp*0.0153; // Last correct
+      //  Float_t digitEn = Calibrate(digitAmp,digitID);
+
+      Float_t etaD=-10., phiD=-10.;
+      fGeom->EtaPhiFromIndex(digitID,etaD,phiD); 
+      //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
+
+      phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
+      
+      Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
+
+      AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
+
+      if(uArray->GetUnitEnergy() == 0.) goodDigit++;
+      uArray->SetUnitTrackID(digitID);
+      
+      Float_t unitEnergy = 0.;
+      Bool_t ok = kFALSE;
+      unitEnergy = uArray->GetUnitEnergy();
+
+      if(unitEnergy==0){
+       if(!fProcId){
+         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+         fProcId = kTRUE;
+       }
+       fRefArray->Add(uArray);
+       fNDigitEmcal++;
+       ok = kTRUE;
+      }
+
+      // Detector flag
+      if(unitEnergy>0.)
+       uArray->SetUnitDetectorFlag(kAll);
+      else uArray->SetUnitDetectorFlag(kEmcal);
+      
+      uArray->SetUnitEnergy(unitEnergy+digitEt);
+
+      uArray->SetUnitCutFlag(kPtHigher);
+
+      // To be modified !!!
+      uArray->SetUnitSignalFlag(kGood);
+       
+      // This is for jet multiplicity
+      uArray->SetUnitClusterID(index);
+       
+      if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
+      
+    } // End loop over cells
+  } // end if !fCluster
+  else { // Keep digit information from clusterization
+
+    // Loop over calo clusters
+    //------------------------------------------------------------------
+
+    // select EMCAL clusters only
+    TRefArray * caloClusters  = new TRefArray();
+    fAOD->GetEMCALClusters(caloClusters);
+
+    // Total number of EMCAL cluster
+    Int_t nclus = caloClusters->GetEntries() ;
+    Int_t beg   = 0;
+    Float_t pos[3] ;
+
+    // Get CaloCells
+    AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
+//(not used ?)    Int_t ncell = cells.GetNumberOfCells() ;
+//(not used ?)    Int_t type = cells.GetType();
+
+    for(Int_t j = beg; j < nclus; j++) { // loop over clusters
+      // Retrieve cluster from aod
+      AliAODCaloCluster *fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
+
+      // Get the cluster info
+//(not used ?)      Float_t energy = fClus->E() ;
+
+      fClus->GetPosition(pos) ;
+      TVector3 vpos(pos[0],pos[1],pos[2]) ;
+
+      Int_t     digMult = fClus->GetNCells() ;
+      UShort_t *digID   = fClus->GetCellsAbsId() ;
+//(not used ?)      Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
+      Int_t     trackIndex = -1;
+      if(fClus->GetNTracksMatched()!=0) 
+       trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();
+//(not used ?)      Int_t     cLabel = fClus->GetLabel(0);
+      // Do double-counted electron correction 
+      if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
+       {
+          // The electron correction go there
+          // Under construction !!!!
+       }  // End of Electron correction
+
+      // Get CaloCells of cluster and fill the unitArray
+      for(Int_t i = 0; i < digMult ; i++) {
+        Int_t    digitID      = digID[i]; // or clus->GetCellNumber(i) ;
+//(not used ?)        Double_t digitAmpFrac = digAmpFrac[i];
+        Float_t  digitAmp     = cells.GetCellAmplitude(digitID) ;
+
+        // Calibration for an energy in GeV
+        Float_t digitEn = digitAmp*0.0153;
+
+        Float_t etaD=-10., phiD=-10.;
+        fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
+        //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
+
+        phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
+
+        Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
+
+       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
+       if(uArray->GetUnitEnergy() == 0.) goodDigit++;
+       uArray->SetUnitTrackID(digitID);
+
+       Float_t unitEnergy = 0.;
+       Bool_t ok = kFALSE;
+       unitEnergy = uArray->GetUnitEnergy();
+
+       if(unitEnergy==0){
+         if(!fProcId){
+           new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+           fProcId = kTRUE;
+         }
+         fRefArray->Add(uArray);
+         fNDigitEmcal++;
+         ok = kTRUE;
+       }
+
+       // Detector flag
+       if(unitEnergy>0.)
+         uArray->SetUnitDetectorFlag(kAll);
+       else uArray->SetUnitDetectorFlag(kEmcal);
+      
+       uArray->SetUnitEnergy(unitEnergy+digitEt);
+
+       uArray->SetUnitCutFlag(kPtHigher);
+
+       // Signal or background
+       // To be modified !!!
+       uArray->SetUnitSignalFlag(kGood);
+
+       // This is for jet multiplicity
+       uArray->SetUnitClusterID(index);
+       
+       if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
+
+      }// End loop over cells
+    } // End loop over clusters
+  } // end else
+
+  fNIn += goodDigit;
+
+  if(fDebug>1)
+    {
+      printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
+      printf("goodDigit : %d\n", goodDigit);
+    }
+}
+
+/*
+//____________________________________________________________________________
+void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
+{
+  // Set calibration parameters:
+  // if calibration database exists, they are read from database,
+  // otherwise, they are taken from digitizer.
+  //
+  // It is a user responsilibity to open CDB before reconstruction,
+  // for example:
+  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
+
+  //Check if calibration is stored in data base
+
+  if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
+    {
+      AliCDBEntry *entry = (AliCDBEntry*)
+        AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+      if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
+    }
+
+  if(!fCalibData)
+    printf("************* Calibration parameters not found in CDB! ****************");
+//    AliFatal("Calibration parameters not found in CDB!");
+
+
+}
+
+//____________________________________________________________________________
+Float_t  AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
+{
+
+  // Convert digitized amplitude into energy.
+  // Calibration parameters are taken from calibration data base for raw data,
+  // or from digitizer parameters for simulated data.
+
+  if(fCalibData){
+
+    if (fGeom==0)
+      printf("************* Did not get geometry from EMCALLoader ***************");
+//      AliFatal("Did not get geometry from EMCALLoader") ;
+
+    Int_t iSupMod = -1;
+    Int_t nModule  = -1;
+    Int_t nIphi   = -1;
+    Int_t nIeta   = -1;
+    Int_t iphi    = -1;
+    Int_t ieta    = -1;
+
+    Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
+    if(!bCell) {
+      // fGeom->PrintGeometry();
+      Error("Calibrate()"," Wrong cell id number : %i", AbsId);
+      assert(0);
+    }
+
+    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
+
+    fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
+    fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
+
+   return -fADCpedestalECA + amp * fADCchannelECA ;
+
+  }
+  else //Return energy with default parameters if calibration is not available
+    return -fADCpedestalECA + amp * fADCchannelECA ;
+
+}
+*/
+
+
diff --git a/JETAN/AliJetAODFillUnitArrayEMCalDigits.h b/JETAN/AliJetAODFillUnitArrayEMCalDigits.h
new file mode 100755 (executable)
index 0000000..eff0287
--- /dev/null
@@ -0,0 +1,76 @@
+#ifndef ALIJETAODFILLUNITARRAYEMCALDIGITS_H
+#define ALIJETAODFILLUNITARRAYEMCALDIGITS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Jet Fill Unit Array 
+// Called by ESD Reader for jet analysis
+// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr)
+//---------------------------------------------------------------------
+
+#include "AliJetFillUnitArray.h"
+
+class AliESDCaloCluster;
+class AliJetReader;
+class AliJetAODReader;
+
+class AliJetAODFillUnitArrayEMCalDigits : public AliJetFillUnitArray
+{
+ public: 
+  AliJetAODFillUnitArrayEMCalDigits();
+  AliJetAODFillUnitArrayEMCalDigits(AliAODEvent *fAOD);
+  virtual ~AliJetAODFillUnitArrayEMCalDigits();
+  
+  // Setter
+  void SetApplyElectronCorrection(Int_t flag = 1)     {fApplyElectronCorrection = flag;}
+  void SetApplyFractionHadronicCorrection(Bool_t val) {fApplyFractionHadronicCorrection = val;}
+  void SetFractionHadronicCorrection(Double_t val)    {fFractionHadronicCorrection = val;}
+  void SetAOD(AliAODEvent* const aod)                 {fAOD = aod;}
+  void SetInitMult(Int_t mult)                        {fNDigitEmcal = mult;}
+  void SetInitMultCut(Int_t multcut)                  {fNDigitEmcalCut = multcut;}
+
+  // Getter
+  Int_t         GetMult()      const {return fNDigitEmcal;}
+  Int_t         GetMultCut()   const {return fNDigitEmcalCut;}
+
+  // For calibration 
+  //  virtual Float_t Calibrate(Int_t amp, Int_t cellId) ;  // Tranforms Amp to energy
+
+  // Other
+  void Exec(Option_t* const option);
+
+ protected:
+  AliAODEvent  *fAOD;                 // ESD
+  Int_t        fNIn;                  // Number of Array filled in UnitArray
+  Int_t        fCluster;              // Use all cells or cells in clusters for jet finding 
+  Int_t        fNCEMCAL;              // Number of clusters in EMCAL
+  Int_t        fNCPHOS;               // Number of clusters in PHOS
+  Int_t        fNCCalo;               // Number of cluster in EMCAL + PHOS calorimeters
+
+  Bool_t       fApplyElectronCorrection;          // Electron correction flag
+  Bool_t       fApplyFractionHadronicCorrection;  // Fraction hadronic correction flag
+  Bool_t       fFractionHadronicCorrection;       // Fraction hadronic correction 
+
+  
+  // geometry info
+  AliESDCaloCluster *fClus;           //! 
+  Int_t fNDigitEmcal;                 //!
+  Int_t fNDigitEmcalCut;              //!
+  //Calibration parameters... to be replaced by database
+//  AliEMCALCalibData *fCalibData;     //! Calibration database if aval
+//  Float_t            fADCchannelECA; // width of one ADC channel for EC section (GeV)
+//  Float_t            fADCpedestalECA;// pedestal of ADC for EC section (GeV)
+
+
+ private:
+  AliJetAODFillUnitArrayEMCalDigits(const AliJetAODFillUnitArrayEMCalDigits &det);
+  AliJetAODFillUnitArrayEMCalDigits &operator=(const AliJetAODFillUnitArrayEMCalDigits &det);
+
+//  void    GetCalibrationParameters(void) ;
+  
+  ClassDef(AliJetAODFillUnitArrayEMCalDigits,1) // Fill Unit Array with tpc and/or emcal information
+};
+
+#endif
diff --git a/JETAN/AliJetAODFillUnitArrayTracks.cxx b/JETAN/AliJetAODFillUnitArrayTracks.cxx
new file mode 100644 (file)
index 0000000..59da52c
--- /dev/null
@@ -0,0 +1,616 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+
+//======================================================================
+// ***July 2006
+// Fill Unit Array class 
+// Class used by AliJetESDReader to fill a UnitArray from the information extracted 
+// from the particle tracks
+// Author: magali.estienne@ires.in2p3.fr
+//======================================================================
+
+
+// --- ROOT system ---
+#include <TVector3.h>
+#include <TProcessID.h>
+
+#include "AliJetUnitArray.h"
+#include "AliJetAODFillUnitArrayTracks.h"
+
+// --- ROOT system ---
+class TSystem;
+class TLorentzVector;
+class TGeoManager;
+
+// --- AliRoot header files ---
+class AliJetFinder;
+
+ClassImp(AliJetAODFillUnitArrayTracks)
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks()
+  : AliJetFillUnitArray(),
+    fNumUnits(0),
+    fHadCorr(0),
+    fApplyMIPCorrection(kTRUE),
+    fAOD(0x0),
+    fGrid0(0x0),
+    fGrid1(0x0),
+    fGrid2(0x0),
+    fGrid3(0x0),
+    fGrid4(0x0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(AliAODEvent* aod)
+  : AliJetFillUnitArray(),
+    fNumUnits(0),
+    fHadCorr(0),
+    fApplyMIPCorrection(kTRUE),
+    fAOD(aod),
+    fGrid0(0x0),
+    fGrid1(0x0),
+    fGrid2(0x0),
+    fGrid3(0x0),
+    fGrid4(0x0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(const AliJetAODFillUnitArrayTracks &det)
+  : AliJetFillUnitArray(det),
+    fNumUnits(det.fNumUnits),
+    fHadCorr(det.fHadCorr),
+    fApplyMIPCorrection(det.fApplyMIPCorrection),
+    fAOD(det.fAOD),
+    fGrid0(det.fGrid0),
+    fGrid1(det.fGrid1),
+    fGrid2(det.fGrid2),
+    fGrid3(det.fGrid3),
+    fGrid4(det.fGrid4)
+{
+  // Copy constructor
+}
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayTracks& AliJetAODFillUnitArrayTracks::operator=(const AliJetAODFillUnitArrayTracks& other)
+{
+  // Assignment
+
+    fNumUnits = other.fNumUnits;
+    fHadCorr = other.fHadCorr;
+    fApplyMIPCorrection = other.fApplyMIPCorrection;
+    fAOD = other.fAOD;
+    fGrid0 = other.fGrid0;
+    fGrid1 = other.fGrid1;
+    fGrid2 = other.fGrid2;
+    fGrid3 = other.fGrid3;
+    fGrid4 = other.fGrid4;
+
+    return (*this);
+}
+
+//____________________________________________________________________________
+void AliJetAODFillUnitArrayTracks::InitParameters()
+{
+  fHadCorr        = 0;     // For hadron correction
+  fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
+
+  fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, 
+                       fPhiMax,fEtaMin,fEtaMax);
+  fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc, 
+                       fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
+
+  fIndex = fTPCGrid->GetIndexObject();
+
+  if(fDebug>20){
+    for(Int_t i=0; i<fNphi+1; i++)
+      for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
+         (*fIndex)(i,j) << endl; }
+  } 
+  if(fDebug>1) printf("\n Parameters initiated ! \n");
+}
+
+//_____________________________________________________________________________
+AliJetAODFillUnitArrayTracks::~AliJetAODFillUnitArrayTracks()
+{
+  // destructor
+}
+
+//_____________________________________________________________________________
+void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/)
+//void AliJetAODFillUnitArrayTracks::Exec(Option_t* option)
+{
+  //
+  // Main method.
+  // Fill the unit array with the charged particle information in AOD
+  //
+
+  fDebug = fReaderHeader->GetDebug();
+  
+  // Set parameters
+  InitParameters();
+
+  // get number of tracks in event (for the loop)
+  Int_t goodTrack = 0;
+  Int_t nt = 0;
+//(not used ?)  Int_t nt2 = 0;
+  Int_t nmax = 0;
+  Float_t pt,eta,phi;
+  //  Int_t sflag = 0;
+  TVector3 p3;
+
+  nt = fAOD->GetNumberOfTracks();
+  if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
+  // 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();  
+  fOpt = fReaderHeader->GetDetector();
+  fDZ  = fReaderHeader->GetDZ();
+  UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
+
+  Int_t nTracksEmcal      = 0;
+  Int_t nTracksEmcalDZ    = 0;
+  Int_t nTracksTpc        = 0;
+  Int_t nTracksTpcOnly    = 0;
+  Int_t nTracksEmcalCut   = 0;
+  Int_t nTracksEmcalDZCut = 0;
+  Int_t nTracksTpcCut     = 0;
+  Int_t nTracksTpcOnlyCut = 0;
+
+//(not used ?)  Int_t nElements  = fTPCGrid->GetNEntries();
+//(not used ?)  Int_t nElements2 = fEMCalGrid->GetNEntries();
+  fGrid = fTPCGrid->GetGridType();
+
+//(not used ?)  Int_t nEntries = fUnitArray->GetEntries();
+
+//(not used ?)  Int_t nRefEnt = fRefArray->GetEntries();
+
+  //loop over tracks
+    nmax = nt;  
+    for (Int_t it = 0; it < nmax; it++) {
+      AliAODTrack *track;
+      track = fAOD->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();
+
+//(not used ?)      Float_t mass = 0.;
+
+      eta = p3.Eta();
+      phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
+      if (status == 0) continue;
+      if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
+
+      if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
+      
+      sflag[goodTrack]=0;
+      if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+      cflag[goodTrack]=0;
+      if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
+
+      if(fGrid==0)
+       {
+         // Only TPC filled from its grid in its total acceptance
+         
+         Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+         Bool_t ok = kFALSE;
+         Bool_t ref = kFALSE;
+
+         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
+         TRefArray *reference = uArray->GetUnitTrackRef();
+         reference->Add(track);
+
+         Float_t unitEnergy = 0.;
+         unitEnergy = uArray->GetUnitEnergy();
+         // nTracksTpcOnly is necessary to count the number of candidate cells
+         // but it doesn't give the real multiplicity -> it will be extracted 
+         // in the jet finder through the reference to tracks
+         if(unitEnergy==0.){
+           nTracksTpcOnly++;
+           ok = kTRUE;
+            ref = kTRUE;
+         }
+
+         // Fill energy in TPC acceptance
+         uArray->SetUnitEnergy(unitEnergy + pt);
+
+         // Pt cut flag
+         if(uArray->GetUnitEnergy()<ptMin){
+           uArray->SetUnitCutFlag(kPtSmaller);
+         }
+         else {
+           uArray->SetUnitCutFlag(kPtHigher);
+           if(ok) nTracksTpcOnlyCut++;
+         }
+
+         // Signal flag
+         if(sflag[goodTrack] == 1) { 
+           uArray->SetUnitSignalFlag(kGood);
+            uArray->SetUnitSignalFlagC(kFALSE,kGood);
+          } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
+
+         if(uArray->GetUnitEnergy()>0 && ref){
+           if(!fProcId) {
+             fProcId = kTRUE;
+             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+           }
+           fRefArray->Add(uArray);
+         }
+       }
+    
+      if(fGrid==1)
+       {
+
+         Int_t nElements = fTPCGrid->GetNEntries();
+         // Fill track information in EMCAL acceptance
+         if((eta >= fEtaMin && eta <= fEtaMax) &&
+            (phi >= fPhiMin && phi <= fPhiMax))// &&
+           {
+             // Include dead-zones
+             if(fDZ)
+               {
+                 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);
+                 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
+                 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
+                 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
+                 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
+                 Int_t n0 = fGrid0->GetNEntries();
+                 Int_t n1 = fGrid1->GetNEntries();
+                 Int_t n2 = fGrid2->GetNEntries();
+                 Int_t n3 = fGrid3->GetNEntries();
+//(not used ?)           Int_t n4 = fGrid4->GetNEntries();
+                 
+                 //              cout << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl;
+                 if(phi >= phimin0 && phi <= phimax0){
+                   Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
+                   AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
+                   TRefArray *reference0 = uArray0->GetUnitTrackRef();
+                   reference0->Add(track);
+                   uArray0->SetUnitTrackID(it);
+
+                   Float_t uEnergy0 = uArray0->GetUnitEnergy();
+                   Bool_t ok0 = kFALSE;
+                   Bool_t ref0 = kFALSE;
+                   if(uEnergy0==0.){
+                     nTracksEmcalDZ++;
+                     ok0 = kTRUE;
+                     ref0 = kTRUE;
+                   }
+                   uArray0->SetUnitEnergy(uEnergy0+pt);
+
+                   if(uArray0->GetUnitEnergy()<ptMin)
+                     uArray0->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray0->SetUnitCutFlag(kPtHigher);
+                     if(ok0) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray0->SetUnitSignalFlag(kGood);
+                      uArray0->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray0->GetUnitEnergy()>0 && ref0){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray0);
+                   }
+                 }
+
+                 if(phi >= phimin1 && phi <= phimax1){
+                   Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
+                   AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
+                   TRefArray *reference1 = uArray1->GetUnitTrackRef();
+                   reference1->Add(track);
+
+                   Float_t uEnergy1 = uArray1->GetUnitEnergy();
+                   Bool_t ok1 = kFALSE;
+                   Bool_t ref1 = kFALSE;
+                   if(uEnergy1==0.){
+                     nTracksEmcalDZ++;
+                     ok1 = kTRUE;
+                     ref1 = kTRUE;
+                   }
+                   uArray1->SetUnitEnergy(uEnergy1+pt);
+
+                   if(uArray1->GetUnitEnergy()<ptMin)
+                     uArray1->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray1->SetUnitCutFlag(kPtHigher);
+                     if(ok1) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray1->SetUnitSignalFlag(kGood);
+                      uArray1->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray1->GetUnitEnergy()>0 && ref1){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray1);
+                   }
+                 }
+
+                 if(phi >= phimin2 && phi <= phimax2){
+                   Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
+                   AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
+                   TRefArray *reference2 = uArray2->GetUnitTrackRef();
+                   reference2->Add(track);
+
+                   Float_t uEnergy2 = uArray2->GetUnitEnergy();
+                   Bool_t ok2 = kFALSE;
+                   Bool_t ref2 = kFALSE;
+                   if(uEnergy2==0.){
+                     nTracksEmcalDZ++;
+                     ok2 = kTRUE;
+                     ref2 = kTRUE;
+                   }
+                   uArray2->SetUnitEnergy(uEnergy2+pt);
+
+                   if(uArray2->GetUnitEnergy()<ptMin)
+                     uArray2->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray2->SetUnitCutFlag(kPtHigher);
+                     if(ok2) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray2->SetUnitSignalFlag(kGood);
+                      uArray2->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray2->GetUnitEnergy()>0 && ref2){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray2);
+                   }
+                 }
+
+                 if(phi >= phimin3 && phi <= phimax3){
+                   Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
+                   AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
+                   TRefArray *reference3 = uArray3->GetUnitTrackRef();
+                   reference3->Add(track);
+
+                   Float_t uEnergy3 = uArray3->GetUnitEnergy();
+                   Bool_t ok3 = kFALSE;
+                   Bool_t ref3 = kFALSE;
+                   if(uEnergy3==0.){
+                     nTracksEmcalDZ++;
+                     ok3 = kTRUE;
+                     ref3 = kTRUE;
+                   }
+                   uArray3->SetUnitEnergy(uEnergy3+pt);
+
+                   if(uArray3->GetUnitEnergy()<ptMin)
+                     uArray3->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray3->SetUnitCutFlag(kPtHigher);
+                     if(ok3) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray3->SetUnitSignalFlag(kGood);
+                      uArray3->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray3->GetUnitEnergy()>0 && ref3){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray3);
+                   }
+                 }
+
+                 if(phi >= phimin4 && phi <= phimax4){
+                   Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
+                   AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
+                   TRefArray *reference4 = uArray4->GetUnitTrackRef();
+                   reference4->Add(track);
+
+                   Float_t uEnergy4 = uArray4->GetUnitEnergy();
+                   Bool_t ok4 = kFALSE;
+                   Bool_t ref4 = kFALSE;
+                   if(uEnergy4==0.){
+                     nTracksEmcalDZ++;
+                     ok4 = kTRUE;
+                     ref4 = kTRUE;
+                   }
+                   uArray4->SetUnitEnergy(uEnergy4+pt);
+
+                   if(uArray4->GetUnitEnergy()<ptMin)
+                     uArray4->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray4->SetUnitCutFlag(kPtHigher);
+                     if(ok4) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray4->SetUnitSignalFlag(kGood);
+                      uArray4->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray4->GetUnitEnergy()>0 && ref4){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray4);
+                   }
+                 }
+               } // end fDZ
+           
+             Int_t towerID = 0;
+             fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
+             if(towerID==-1) continue;
+             
+             AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
+             TRefArray *reference = uArray->GetUnitTrackRef();
+             reference->Add(track);
+
+             Float_t unitEnergy = uArray->GetUnitEnergy(); 
+             Bool_t ok = kFALSE;
+             Bool_t ref = kFALSE;
+             if(unitEnergy==0.){
+               nTracksEmcal++;
+               ok=kTRUE;
+               ref=kTRUE;
+             }
+
+             // Do Hadron Correction
+             // This is under construction !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+             // Parametrization to be added
+             if (fApplyMIPCorrection != 0) 
+               { 
+//               Float_t   hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
+//               unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
+               } //end Hadron Correction loop
+
+             uArray->SetUnitEnergy(unitEnergy + pt);
+
+             // Put a pt cut flag
+             if(uArray->GetUnitEnergy()<ptMin){
+               uArray->SetUnitCutFlag(kPtSmaller);
+             }
+             else {
+               uArray->SetUnitCutFlag(kPtHigher);
+               if(ok) nTracksEmcalCut++;
+             }
+
+             // Signal flag
+             if(sflag[goodTrack] == 1) { 
+               uArray->SetUnitSignalFlag(kGood);
+                uArray->SetUnitSignalFlagC(kFALSE,kGood);
+              } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
+
+             if(uArray->GetUnitEnergy()>0 && ref){
+               if(!fProcId){
+                 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+                 fProcId = kTRUE;
+               }
+               fRefArray->Add(uArray);
+             }
+           } // end loop on EMCal acceptance cut + tracks quality
+         else{ 
+           // Outside EMCal acceptance
+           //  Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
+           Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+
+           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
+           TRefArray *reference = uArray->GetUnitTrackRef();
+           reference->Add(track);
+
+           Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
+           Bool_t okout = kFALSE;
+           Bool_t refout = kFALSE;
+           if(unitEnergy2==0.){
+             nTracksTpc++;
+             okout=kTRUE;
+             refout=kTRUE;
+           }
+           // Fill energy outside emcal acceptance
+           uArray->SetUnitEnergy(unitEnergy2 + pt);
+         
+           // Pt cut flag
+           if(uArray->GetUnitEnergy()<ptMin){
+             uArray->SetUnitCutFlag(kPtSmaller);
+           }
+           else {
+             uArray->SetUnitCutFlag(kPtHigher);
+             if(okout) nTracksTpcCut++;
+           }
+           // Signal flag
+           if(sflag[goodTrack] == 1) {
+             uArray->SetUnitSignalFlag(kGood);
+              uArray->SetUnitSignalFlagC(kFALSE,kGood);
+            } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
+
+           if(uArray->GetUnitEnergy()>0 && refout){
+             if(!fProcId){
+               new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+               fProcId = kTRUE;
+             }
+             fRefArray->Add(uArray);
+           }
+         }
+       } // end fGrid==1
+
+      goodTrack++;
+
+    } // end loop on entries (tpc tracks)
+
+
+  if(fDebug>0)
+    {
+      cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
+      cout << "goodTracks: " << goodTrack << endl;
+    }
+
+  delete[] sflag;
+  delete[] cflag;
+
+  //    } // end loop on entries (tpc tracks)
+      
+  if(fGrid==0) {
+    fNTracks = nTracksTpcOnly;
+    fNTracksCut = nTracksTpcOnlyCut;
+    if(fDebug>10){
+      cout << "fNTracks : " << fNTracks << endl;
+      cout << "fNTracksCut : " << fNTracksCut << endl;
+    }
+  }
+  if(fGrid==1) {
+    fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
+    fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
+    if(fDebug>10){
+      cout << "fNTracks : " << fNTracks << endl;
+      cout << "fNTracksCut : " << fNTracksCut << endl;
+    }
+  }  
+  
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/JETAN/AliJetAODFillUnitArrayTracks.h b/JETAN/AliJetAODFillUnitArrayTracks.h
new file mode 100644 (file)
index 0000000..67dd8fa
--- /dev/null
@@ -0,0 +1,64 @@
+#ifndef ALIJETAODFILLUNITARRAYTRACKS_H
+#define ALIJETAODFILLUNITARRAYTRACKS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Jet Fill Unit Array 
+// Called by ESD Reader for jet analysis
+// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr)
+//---------------------------------------------------------------------
+
+#include "AliJetFillUnitArray.h"
+
+class AliJetReader;
+class AliJetAODReader;
+class AliEMCALGeometry;
+
+class AliJetAODFillUnitArrayTracks : public AliJetFillUnitArray
+{
+ public: 
+  AliJetAODFillUnitArrayTracks();
+  AliJetAODFillUnitArrayTracks(AliAODEvent *fAOD);
+  virtual ~AliJetAODFillUnitArrayTracks();
+  
+  // Setter
+  void SetHadCorrector(AliJetHadronCorrection* const corr) {fHadCorr = corr;}
+  void SetApplyMIPCorrection(Bool_t const val)             {fApplyMIPCorrection = val;}
+  void SetAOD(AliAODEvent* const aod)                      {fAOD = aod;}
+  void SetGrid0(AliJetGrid* const grid0)                   {fGrid0 = grid0;}
+  void SetGrid1(AliJetGrid* const grid1)                   {fGrid1 = grid1;}
+  void SetGrid2(AliJetGrid* const grid2)                   {fGrid2 = grid2;}
+  void SetGrid3(AliJetGrid* const grid3)                   {fGrid3 = grid3;}
+  void SetGrid4(AliJetGrid* const grid4)                   {fGrid4 = grid4;}
+
+  // Getter
+  Int_t GetHadCorrection()  const {return fApplyMIPCorrection;}
+  Int_t GetMult()           const {return fNTracks;}
+  Int_t GetMultCut()        const {return fNTracksCut;}
+
+  // Other
+  void Exec(Option_t* const option);
+
+ protected:
+  Int_t                     fNumUnits;           // Number of units in the unit object array (same as num towers in EMCAL)
+  AliJetHadronCorrection   *fHadCorr;            // Pointer to Hadron Correction Object
+  Bool_t                    fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection
+
+  AliAODEvent              *fAOD;                // ESD
+  AliJetGrid               *fGrid0;              // Grid used for dead zones definition
+  AliJetGrid               *fGrid1;              // Grid used for dead zones definition
+  AliJetGrid               *fGrid2;              // Grid used for dead zones definition
+  AliJetGrid               *fGrid3;              // Grid used for dead zones definition
+  AliJetGrid               *fGrid4;              // Grid used for dead zones definition
+
+ private:
+  AliJetAODFillUnitArrayTracks(const AliJetAODFillUnitArrayTracks &det);
+  AliJetAODFillUnitArrayTracks &operator=(const AliJetAODFillUnitArrayTracks &det);
+  void InitParameters();
+
+  ClassDef(AliJetAODFillUnitArrayTracks,1) // Fill Unit Array with tpc and/or emcal information
+};
+
+#endif
index 89c7d62..b62ee7f 100644 (file)
 #include <TLorentzVector.h>
 #include <TVector3.h>
 #include <TChain.h>
+#include <TFile.h>
+#include <TTask.h>
+#include <TGeoManager.h>
 
 #include "AliJetAODReader.h"
 #include "AliJetAODReaderHeader.h"
 #include "AliAODEvent.h"
 #include "AliAODTrack.h"
+#include "AliJetDummyGeo.h"
+#include "AliJetAODFillUnitArrayTracks.h"
+#include "AliJetAODFillUnitArrayEMCalDigits.h"
+#include "AliJetHadronCorrection.h"
+#include "AliJetUnitArray.h"
 
 ClassImp(AliJetAODReader)
 
 AliJetAODReader::AliJetAODReader():
     AliJetReader(),
     fChain(0x0),
+    fTree(0x0),
     fAOD(0x0),
     fRef(new TRefArray),
     fDebug(0),
-    fOpt(0)
+    fOpt(0),
+    fGeom(0),
+    fHadCorr(0x0),
+    fTpcGrid(0x0),
+    fEmcalGrid(0x0),
+    fGrid0(0),
+    fGrid1(0),
+    fGrid2(0),
+    fGrid3(0),
+    fGrid4(0),
+    fPtCut(0),
+    fHCorrection(0),
+    fApplyElectronCorrection(kFALSE),
+    fEFlag(kFALSE),
+    fApplyMIPCorrection(kTRUE),
+    fApplyFractionHadronicCorrection(kFALSE),
+    fFractionHadronicCorrection(0.3),
+    fNumUnits(0),
+    fMass(0),
+    fSign(0),
+    fNIn(0),
+    fDZ(0),
+    fNeta(0),
+    fNphi(0),
+    fArrayInitialised(0),
+    fRefArray(0x0),
+    fProcId(kFALSE)
+
 {
   // Constructor    
 }
@@ -55,6 +91,17 @@ AliJetAODReader::~AliJetAODReader()
     delete fChain;
     delete fAOD;
     delete fRef;
+    delete fTpcGrid;
+    delete fEmcalGrid;
+    if(fDZ)
+      {
+        delete fGrid0;
+        delete fGrid1;
+        delete fGrid2;
+        delete fGrid3;
+        delete fGrid4;
+      }
+
 }
 
 //____________________________________________________________________________
@@ -88,13 +135,12 @@ void AliJetAODReader::OpenInputFiles()
   
   gSystem->FreeDirectory(dir);
   
-
   fAOD = 0;
   fChain->SetBranchAddress("AOD",&fAOD);
   
   int nMax = fChain->GetEntries(); 
 
-  if(fDebug>1)printf("\n AliJetAODReader::OpenInputFiles Total number of events in chain= %d \n",nMax);
+  printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax);
   
   // set number of events in header
   if (fReaderHeader->GetLastEvent() == -1)
@@ -113,7 +159,7 @@ void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) {
      fChain = (TChain*)      tree;
      
      Int_t nMax = fChain->GetEntries(); 
-     if(fDebug>1)printf("\n AliJetAODReader::ConnectTree Total number of events in chain= %5d \n", nMax);
+     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);
@@ -138,7 +184,7 @@ Bool_t AliJetAODReader::FillMomentumArray()
   
   // get number of tracks in event (for the loop)
   Int_t nt = fAOD->GetNTracks();
-  if(fDebug>2)printf(" AliJetAODReader::FillMomentumArray AOD tracks: %5d \t", nt);
+  printf("AOD tracks: %5d \t", nt);
   
   // temporary storage of signal and pt cut flag
   Int_t* sflag  = new Int_t[nt];
@@ -173,7 +219,7 @@ Bool_t AliJetAODReader::FillMomentumArray()
     aodTrack++;
     fRef->Add(track);
   }
-  if(fDebug>2)printf("Used AOD tracks: %5d \n", aodTrack);
+  printf("Used AOD tracks: %5d \n", aodTrack);
   // set the signal flags
   fSignalFlag.Set(aodTrack,sflag);
   fCutFlag.Set(aodTrack,cflag);
@@ -183,3 +229,356 @@ Bool_t AliJetAODReader::FillMomentumArray()
   
   return kTRUE;
 }
+
+//__________________________________________________________
+void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
+{
+  //
+  // Set flag to apply MIP correction fApplyMIPCorrection
+  // - exclusive with fApplyFractionHadronicCorrection
+  //
+
+  fApplyMIPCorrection = val;
+  if(fApplyMIPCorrection == kTRUE)
+    {
+      SetApplyFractionHadronicCorrection(kFALSE);
+      printf("Enabling MIP Correction \n");
+    }
+  else
+    {
+      printf("Disabling MIP Correction \n");
+    }
+}
+
+//__________________________________________________________
+void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
+{
+  //
+  // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
+  // - exclusive with fApplyMIPCorrection
+  //
+
+  fApplyFractionHadronicCorrection = val;
+  if(fApplyFractionHadronicCorrection == kTRUE)
+    {
+      SetApplyMIPCorrection(kFALSE);
+      printf("Enabling Fraction Hadronic Correction \n");
+    }
+  else
+    {
+      printf("Disabling Fraction Hadronic Correction \n");
+    }
+}
+
+//__________________________________________________________
+void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
+{
+  //
+  // Set value to fFractionHadronicCorrection (default is 0.3)
+  // apply EMC hadronic correction fApplyFractionHadronicCorrection
+  // - exclusive with fApplyMIPCorrection
+  //
+
+  fFractionHadronicCorrection = val;
+  if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
+    {
+      SetApplyFractionHadronicCorrection(kTRUE);
+      printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
+    }
+  else
+    {
+      SetApplyFractionHadronicCorrection(kFALSE);
+    }
+}
+
+//____________________________________________________________________________
+void AliJetAODReader::CreateTasks(TChain* tree)
+{
+  fDebug = fReaderHeader->GetDebug();
+  fDZ = fReaderHeader->GetDZ();
+  fTree = tree;
+
+  // Init EMCAL geometry and create UnitArray object
+  SetEMCALGeometry();
+  //  cout << "In create task" << endl;
+  InitParameters();
+  InitUnitArray();
+
+  fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
+  fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
+  fFillUAFromTracks->SetReaderHeader(fReaderHeader);
+  fFillUAFromTracks->SetGeom(fGeom);
+  fFillUAFromTracks->SetTPCGrid(fTpcGrid);
+  fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
+
+  if(fDZ)
+    {
+      fFillUAFromTracks->SetGrid0(fGrid0);
+      fFillUAFromTracks->SetGrid1(fGrid1);
+      fFillUAFromTracks->SetGrid2(fGrid2);
+      fFillUAFromTracks->SetGrid3(fGrid3);
+      fFillUAFromTracks->SetGrid4(fGrid4);
+    }
+  fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
+  fFillUAFromTracks->SetHadCorrector(fHadCorr);
+  fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
+  fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
+  fFillUAFromEMCalDigits->SetGeom(fGeom);
+  fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
+  fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
+  fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
+  fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
+  fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
+
+  fFillUnitArray->Add(fFillUAFromTracks);
+  fFillUnitArray->Add(fFillUAFromEMCalDigits);
+  fFillUAFromTracks->SetActive(kFALSE);
+  fFillUAFromEMCalDigits->SetActive(kFALSE);
+
+  cout << "Tasks instantiated at that stage ! " << endl;
+  cout << "You can loop over events now ! " << endl;
+
+}
+
+//____________________________________________________________________________
+Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
+{
+  //
+  // Main function
+  // Fill the reader part
+  //
+
+  fProcId = procid;
+  fRefArray = refArray;
+//(not used ?)  Int_t nEntRef = fRefArray->GetEntries();
+//(not used ?)  Int_t nEntUnit = fUnitArray->GetEntries();
+
+  // clear momentum array
+  ClearArray();
+
+  fDebug = fReaderHeader->GetDebug();
+  fOpt = fReaderHeader->GetDetector();
+
+  if(!fAOD) {
+    return kFALSE;
+  }
+
+  // TPC only or Digits+TPC or Clusters+TPC
+  if(fOpt%2==!0 && fOpt!=0){
+    fFillUAFromTracks->SetAOD(fAOD);
+    fFillUAFromTracks->SetActive(kTRUE);
+    fFillUAFromTracks->SetUnitArray(fUnitArray);
+    fFillUAFromTracks->SetRefArray(fRefArray);
+    fFillUAFromTracks->SetProcId(fProcId);
+    //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
+    fFillUAFromTracks->Exec("tpc");
+    if(fOpt==1){
+      fNumCandidate = fFillUAFromTracks->GetMult();
+      fNumCandidateCut = fFillUAFromTracks->GetMultCut();
+    }
+  }
+
+  // Digits only or Digits+TPC
+  if(fOpt>=2 && fOpt<=3){
+    fFillUAFromEMCalDigits->SetAOD(fAOD);
+    fFillUAFromEMCalDigits->SetActive(kTRUE);
+    fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
+    fFillUAFromEMCalDigits->SetRefArray(fRefArray);
+    fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
+    fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
+    fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
+    fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
+    fNumCandidate = fFillUAFromEMCalDigits->GetMult();
+    fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
+  }
+
+  //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
+
+  return kTRUE;
+}
+
+//____________________________________________________________________________
+Bool_t AliJetAODReader::SetEMCALGeometry()
+{
+  // 
+  // Set the EMCal Geometry
+  //
+  
+  if (!fTree->GetFile()) 
+    return kFALSE;
+
+  TString geomFile(fTree->GetFile()->GetName());
+  geomFile.ReplaceAll("AliESDs", "geometry");
+  
+  // temporary workaround for PROOF bug #18505
+  geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
+  if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
+
+  // Define EMCAL geometry to be able to read ESDs
+  fGeom = AliJetDummyGeo::GetInstance();
+  if (fGeom == 0)
+    fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
+  
+  // To be setted to run some AliEMCALGeometry functions
+  TGeoManager::Import(geomFile);
+  fGeom->GetTransformationForSM();  
+  printf("\n EMCal Geometry set ! \n");
+  
+  return kTRUE;
+  
+}
+
+//____________________________________________________________________________
+void AliJetAODReader::InitParameters()
+{
+  // Initialise parameters
+  fOpt = fReaderHeader->GetDetector();
+  //  fHCorrection    = 0;          // For hadron correction
+  fHadCorr        = 0;          // For hadron correction
+  if(fEFlag==kFALSE){
+    if(fOpt==0 || fOpt==1)
+      fECorrection    = 0;        // For electron correction
+    else fECorrection = 1;        // For electron correction
+  }
+  fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
+  if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
+}
+
+//____________________________________________________________________________
+void AliJetAODReader::InitUnitArray()
+{
+  //Initialises unit arrays
+  Int_t nElements = fTpcGrid->GetNEntries();
+  Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
+  if(fArrayInitialised) fUnitArray->Delete();
+
+  if(fTpcGrid->GetGridType()==0)
+    { // Fill the following quantities :
+      // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi,
+      // detector flag, in/out jet, pt cut, mass, cluster ID)
+      for(Int_t nBin = 1; nBin < nElements+1; nBin++)
+        {
+          //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
+          fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
+          phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
+          Deta = fTpcGrid->GetDeta();
+          Dphi = fTpcGrid->GetDphi();
+          new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+        }
+    }
+
+  if(fTpcGrid->GetGridType()==1)
+    {
+      Int_t nGaps = 0;
+      Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
+
+      if(fDZ)
+        {
+          // Define a grid of cell for the gaps between SM
+          Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
+          Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
+          fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
+          fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
+          fGrid0->SetGridType(0);
+          fGrid0->SetMatrixIndexes();
+          fGrid0->SetIndexIJ();
+          n0 = fGrid0->GetNEntries();
+          fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
+          fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
+          fGrid1->SetGridType(0);
+          fGrid1->SetMatrixIndexes();
+          fGrid1->SetIndexIJ();
+          n1 = fGrid1->GetNEntries();
+          fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
+          fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
+          fGrid2->SetGridType(0);
+          fGrid2->SetMatrixIndexes();
+          fGrid2->SetIndexIJ();
+          n2 = fGrid2->GetNEntries();
+          fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
+          fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
+          fGrid3->SetGridType(0);
+          fGrid3->SetMatrixIndexes();
+          fGrid3->SetIndexIJ();
+          n3 = fGrid3->GetNEntries();
+          fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
+          fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
+          fGrid4->SetGridType(0);
+          fGrid4->SetMatrixIndexes();
+          fGrid4->SetIndexIJ();
+          n4 = fGrid4->GetNEntries();
+
+          nGaps = n0+n1+n2+n3+n4;
+
+        }
+
+      for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++) 
+       {
+         if(nBin<fNumUnits)
+           {
+             fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry 
+             // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
+             phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
+             Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
+             Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
+             new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+           } 
+         else {
+           if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
+             fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
+             phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
+             Deta = fTpcGrid->GetDeta();
+             Dphi = fTpcGrid->GetDphi();
+             new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+           }
+           else {
+             if(fDZ) {
+               if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
+                 if(nBin<fNumUnits+nElements+n0)
+                   {
+                     Float_t phi = eta = 0.;
+                     fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
+                     Deta = fGrid0->GetDeta(); 
+                     Dphi = fGrid0->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                   }
+                 else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
+                   {
+                     Float_t phi = eta = 0.;
+                     fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
+                     Deta = fGrid1->GetDeta(); 
+                     Dphi = fGrid1->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                   }
+                 else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
+                   {
+                     Float_t phi = eta = 0.;
+                     fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
+                     Deta = fGrid2->GetDeta(); 
+                     Dphi = fGrid2->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                   }
+                 else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
+                   {
+                     Float_t phi = eta = 0.;
+                     fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
+                     Deta = fGrid3->GetDeta(); 
+                     Dphi = fGrid3->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                   }
+                 else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
+                   {
+                     Float_t phi = eta = 0.;
+                     fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
+                     Deta = fGrid4->GetDeta(); 
+                     Dphi = fGrid4->GetDphi(); 
+                     new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
+                   }
+               }
+             } // end if(fDZ)
+           } // end else 2
+         } // end else 1
+       } // end loop on nBin
+    } // end grid type == 1
+  fArrayInitialised = 1;
+}
index 889f15a..d11bbb2 100644 (file)
 //---------------------------------------------------------------------
 
 #include "AliJetReader.h"
+#include "AliJetUnitArray.h"
+#include "AliJetGrid.h"
+class AliJetUnitArray;
+class AliJetDummyGeo;
+class AliJetHadronCorrection;
 class AliJetAODReaderHeader;
 class AliJetReaderHeader;
 class AliAODEvent;
@@ -27,17 +32,66 @@ class AliJetAODReader : public AliJetReader
   Bool_t FillMomentumArray(); 
   void   OpenInputFiles();
   void   ConnectTree(TTree* tree, TObject* data);
-  void   SetInputEvent(TObject* /*esd*/, TObject* aod, TObject* /*mc*/) {fAOD = (AliAODEvent*) aod;}
+  void   InitUnitArray();
+  void   CreateTasks(TChain* tree);
+  Bool_t ExecTasks(Bool_t procid, TRefArray* refArray);
+
+  void SetInputEvent(TObject* /*esd*/, TObject* aod, TObject* /*mc*/) {fAOD = (AliAODEvent*) aod;}
+  void SetTPCGrid(AliJetGrid *grid)   {fTpcGrid = grid;}
+  void SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;}
+  // Correction of hadronic energy
+  void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;}
+  void SetHadronCorrector(AliJetHadronCorrection* corr) {fHadCorr = corr;}
+  void SetApplyElectronCorrection(Int_t flag = 1) {fECorrection = flag; fEFlag=kTRUE;}
+  void SetApplyMIPCorrection(Bool_t val);
+  void SetApplyFractionHadronicCorrection(Bool_t val);
+  void SetFractionHadronicCorrection(Double_t val);
+
+
  private:
+  Bool_t SetEMCALGeometry();
+  void InitParameters();
   AliJetAODReader(const AliJetAODReader &det);
   AliJetAODReader &operator=(const AliJetAODReader &det);
 
  private:
   TChain                     *fChain;  //! chain for reconstructed tracks
+  TChain                     *fTree;   //! tree to get the EMCal geometry
   AliAODEvent                *fAOD;    //! pointer to aod
-  TRefArray                  *fRef;    //! pointer to array of references to tracks
+  TRefArray                  *fRef;    // pointer to array of references to tracks
   Int_t                       fDebug;  // Debug option
   Int_t                       fOpt;    // Detector to be used for jet reconstruction
+  AliJetDummyGeo             *fGeom;             //! EMCAL Geometry
+
+  AliJetHadronCorrection     *fHadCorr;          //! Pointer to Hadron Correction Object
+  AliJetGrid                 *fTpcGrid;          //! Pointer to grid object
+  AliJetGrid                 *fEmcalGrid;        //! Pointer to grid object
+  AliJetGrid                 *fGrid0;            // Pointer to grid object
+  AliJetGrid                 *fGrid1;            // Pointer to grid object
+  AliJetGrid                 *fGrid2;            // Pointer to grid object
+  AliJetGrid                 *fGrid3;            // Pointer to grid object
+  AliJetGrid                 *fGrid4;            // Pointer to grid object
+  Float_t                     fPtCut;            // Pt cut for tracks to minimise background contribution
+  Int_t                       fHCorrection;      // Hadron correction flag
+  Int_t                       fApplyElectronCorrection;      // Electron correction flag
+  Bool_t                      fEFlag;            // if (fEFlag == kFALSE) => fECorrection automatically setted
+  Bool_t                      fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection
+  Bool_t                      fApplyFractionHadronicCorrection; // Another type of charged particle energy deposition in EMC
+  Double_t                    fFractionHadronicCorrection; // Fraction of momentum of the TPC track to be subtracted from EMC tower
+  Int_t                       fNumUnits;         // Number of units in the unit object array
+                                                 // (same as num towers in EMCAL)
+  Float_t                     fMass;             // Particle mass
+  Int_t                       fSign;             // Particle sign
+  Int_t                       fNIn;              // Number of Array filled in UnitArray
+  Bool_t                      fDZ;               // Use or not dead zones
+  Int_t                       fNeta;             // Number of bins in eta of tpc grid
+  Int_t                       fNphi;             // Number of bins in phi of tpc grid
+  Bool_t                      fArrayInitialised; // To check that array of units is initialised
+
+  TRefArray                  *fRefArray;         // array of digit position and energy
+  Bool_t                      fProcId;           // Bool_t for TProcessID synchronization
+
+
   ClassDef(AliJetAODReader,1)
 };
  
diff --git a/JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx b/JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx
new file mode 100755 (executable)
index 0000000..199d64d
--- /dev/null
@@ -0,0 +1,415 @@
+
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+//---------------------------------------------------------------------
+// Fill Unit Array 
+// Called by ESD reader for jet analysis
+// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
+//---------------------------------------------------------------------
+
+// --- ROOT system ---
+#include <TLorentzVector.h>
+#include <TMath.h>
+
+// --- AliRoot header files ---
+#include "AliESDCaloCluster.h"
+#include "AliESDCaloCells.h"
+#include "AliJetUnitArray.h"
+#include "AliJetESDFillUnitArrayEMCalDigits.h"
+// #include "AliEMCALCalibData.h"
+// #include "AliCDBManager.h"
+// // class AliCDBStorage;
+// #include "AliCDBEntry.h"
+
+// --- ROOT system ---
+class TSystem;
+class TGeoManager;
+class TArrayS;
+
+// --- AliRoot header files ---
+class AliJetFinder;
+class AliJetReader;
+class AliJetESDReader;
+// class AliEMCALCalibData;
+// class AliCDBManager;
+// class AliCDBStorage;
+// class AliCDBEntry;
+
+ClassImp(AliJetESDFillUnitArrayEMCalDigits)
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits()
+  : AliJetFillUnitArray(),
+    fESD(0),
+    fNIn(0),
+    fOpt(0),
+    fCluster(0),
+    fDebug(0),
+    fNCEMCAL(0),
+    fNCPHOS(0),
+    fNCCalo(0),
+    fApplyElectronCorrection(kFALSE),
+    fApplyFractionHadronicCorrection(kFALSE),
+    fFractionHadronicCorrection(0.3),
+    fClus(0x0),
+    fNDigitEmcal(0),
+    fNDigitEmcalCut(0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *esd)
+  : AliJetFillUnitArray(),
+    fESD(esd),
+    fNIn(0),
+    fOpt(0),
+    fCluster(0),
+    fDebug(0),
+    fNCEMCAL(0),
+    fNCPHOS(0),
+    fNCCalo(0),
+    fApplyElectronCorrection(kFALSE),
+    fApplyFractionHadronicCorrection(kFALSE),
+    fFractionHadronicCorrection(0.3),
+    fClus(0x0),
+    fNDigitEmcal(0),
+    fNDigitEmcalCut(0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det)
+  : AliJetFillUnitArray(det),
+    fESD(det.fESD),
+    fNIn(det.fNIn),
+    fOpt(det.fOpt),
+    fCluster(det.fCluster),
+    fDebug(det.fDebug),
+    fNCEMCAL(det.fNCEMCAL),
+    fNCPHOS(det.fNCPHOS),
+    fNCCalo(det.fNCCalo),
+    fApplyElectronCorrection(det.fApplyElectronCorrection),
+    fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
+    fFractionHadronicCorrection(det.fFractionHadronicCorrection),
+    fClus(det.fClus),
+    fNDigitEmcal(det.fNDigitEmcal),
+    fNDigitEmcalCut(det.fNDigitEmcalCut)
+{
+  // Copy constructor
+}
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayEMCalDigits& AliJetESDFillUnitArrayEMCalDigits::operator=(const AliJetESDFillUnitArrayEMCalDigits& other)
+{
+  // Assignment
+
+  fESD = other.fESD;
+  fNIn = other.fNIn;
+  fOpt = other.fOpt;
+  fCluster = other.fCluster;
+  fDebug = other.fDebug;
+  fNCEMCAL = other.fNCEMCAL;
+  fNCPHOS = other.fNCPHOS;
+  fNCCalo = other.fNCCalo;
+  fApplyElectronCorrection = other.fApplyElectronCorrection;
+  fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection;
+  fFractionHadronicCorrection = other.fFractionHadronicCorrection;
+  fClus = other.fClus;
+  fNDigitEmcal = other.fNDigitEmcal;
+  fNDigitEmcalCut = other.fNDigitEmcalCut;
+
+  return (*this);
+
+}
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits()
+{
+  // destructor
+}
+
+//_____________________________________________________________________________
+void AliJetESDFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
+{
+  //
+  // Main method.
+  // Fill the unit array with the neutral particle information from the EMCal cells in ESD
+  //
+
+  fDebug = fReaderHeader->GetDebug();
+  fOpt = fReaderHeader->GetDetector();
+  fCluster = fReaderHeader->GetCluster();
+
+  // Init parameters
+  InitParameters();
+
+//(not used ?)  Int_t   nDigitTot      = 0;
+  Int_t   goodDigit      = 0;
+//(not used ?)  Int_t   beg            = 0;
+//(not used ?)  Int_t   end            = 0; 
+  Int_t   index          = 0;
+//(not used ?)  Int_t   count          = 0;
+
+//(not used ?)  Int_t nRefEnt = fRefArray->GetEntries();
+
+  if(!fCluster) { // Keep all digit information
+    // Loop over all cell information
+    //------------------------------------------------------------------
+    AliESDCaloCells &cells= *(fESD->GetEMCALCells());
+    Int_t ncell = cells.GetNumberOfCells() ;
+//(not used ?)    Int_t type = cells.GetType();
+
+    for (Int_t icell=  0; icell <  ncell; icell++) {
+      Int_t   digitID   = cells.GetCellNumber(icell);
+      Float_t digitAmp  = cells.GetAmplitude(icell);
+//(not used ?)      Float_t digitTime = cells.GetTime(icell);
+      Float_t digitEn   = digitAmp*0.0153; // Last correct
+      //  Float_t digitEn = Calibrate(digitAmp,digitID);
+
+      Float_t etaD=-10., phiD=-10.;
+      fGeom->EtaPhiFromIndex(digitID,etaD,phiD); 
+      //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
+
+      phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
+      
+      Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
+
+      AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
+
+      if(uArray->GetUnitEnergy() == 0.) goodDigit++;
+      uArray->SetUnitTrackID(digitID);
+
+      Float_t unitEnergy = 0.;
+      Bool_t ok = kFALSE;
+      unitEnergy = uArray->GetUnitEnergy();
+
+      if(unitEnergy==0){
+       if(!fProcId){
+         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+         fProcId = kTRUE;
+       }
+       fRefArray->Add(uArray);
+       fNDigitEmcal++;
+       ok = kTRUE;
+      }
+
+      // Detector flag
+      if(unitEnergy>0.)
+       uArray->SetUnitDetectorFlag(kAll);
+      else uArray->SetUnitDetectorFlag(kEmcal);
+      
+      uArray->SetUnitEnergy(unitEnergy+digitEt);
+
+      uArray->SetUnitCutFlag(kPtHigher);
+
+      // To be modified !!!
+      uArray->SetUnitSignalFlag(kGood);
+
+      // This is for jet multiplicity
+      uArray->SetUnitClusterID(index);
+       
+      if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
+
+    } // End loop over cells
+  } // end if !fCluster
+  else { // Keep digit information from clusterization
+
+    // Loop over calo clusters
+    //------------------------------------------------------------------
+
+    //select EMCAL clusters only
+    TRefArray * caloClusters  = new TRefArray();
+    fESD->GetEMCALClusters(caloClusters);
+
+    // Total number of EMCAL cluster
+    Int_t nclus = caloClusters->GetEntries() ;
+    Int_t beg   = 0;
+    Float_t pos[3] ;
+
+    // Get reconstructed vertex position
+    Double_t vertexPosition[3] ;
+    fESD->GetVertex()->GetXYZ(vertexPosition) ;
+
+    // Get CaloCells
+    AliESDCaloCells &cells= *(fESD->GetEMCALCells());
+//(not used ?)    Int_t ncell = cells.GetNumberOfCells() ;
+//(not used ?)    Int_t type = cells.GetType();
+
+    for(Int_t j = beg; j < nclus; j++) { // loop over clusters
+      // Retrieve cluster from esd
+      AliESDCaloCluster *fClus = (AliESDCaloCluster *) caloClusters->At(j) ;
+
+      // Get the cluster info
+//(not used ?)      Float_t energy         = fClus->E() ;
+//(not used ?)      Int_t   iprim          = fClus->GetLabel();
+
+      fClus->GetPosition(pos) ;
+      TVector3 vpos(pos[0],pos[1],pos[2]) ;
+      TLorentzVector p ;
+      fClus->GetMomentum(p,vertexPosition);
+
+      Int_t     digMult = fClus->GetNCells() ;
+      UShort_t *digID   = fClus->GetCellsAbsId() ;
+//(not used ?)      Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
+      Int_t     trackIndex = fClus->GetTrackMatched();
+
+      // Do double-counted electron correction 
+      if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
+       {
+          // The electron correction go there
+          // Under construction !!!!
+       }  // End of Electron correction
+
+      // Get CaloCells of cluster and fill the unitArray
+      for(Int_t i = 0; i < digMult ; i++){
+        Int_t    digitID   = digID[i]; // or clus->GetCellNumber(i) ;
+//(not used ?)        Double_t digitAmpFrac = digAmpFrac[i];
+        Float_t  digitAmp     = cells.GetCellAmplitude(digitID) ;
+//(not used ?)        Float_t  digitTime    = cells.GetCellTime(digitID);
+
+        // Calibration for an energy in GeV
+        Float_t digitEn = digitAmp*0.0153;
+
+        Float_t etaD=-10., phiD=-10.;
+        fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
+        //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
+
+        phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
+
+        Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
+
+       cout << "Digit " << i << ", eta: " << etaD << ", phi: " << phiD << endl;
+
+       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
+       if(uArray->GetUnitEnergy() == 0.) goodDigit++;
+       uArray->SetUnitTrackID(digitID);
+
+       Float_t unitEnergy = 0.;
+       Bool_t ok = kFALSE;
+       unitEnergy = uArray->GetUnitEnergy();
+
+       if(unitEnergy==0){
+         if(!fProcId){
+           new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+           fProcId = kTRUE;
+         }
+         fRefArray->Add(uArray);
+         fNDigitEmcal++;
+         ok = kTRUE;
+       }
+
+       // Detector flag
+       if(unitEnergy>0.)
+         uArray->SetUnitDetectorFlag(kAll);
+       else uArray->SetUnitDetectorFlag(kEmcal);
+      
+       uArray->SetUnitEnergy(unitEnergy+digitEt);
+
+       uArray->SetUnitCutFlag(kPtHigher);
+
+       // To be modified !!!
+       uArray->SetUnitSignalFlag(kGood);
+
+       // This is for jet multiplicity
+       uArray->SetUnitClusterID(index);
+       
+       if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
+      
+      } // End loop over cells
+    } // End loop over clusters
+  } // end else
+
+  fNIn += goodDigit;
+
+  if(fDebug>1)
+    {
+      printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
+      printf("goodDigit : %d\n", goodDigit);
+    }
+}
+
+/*
+//____________________________________________________________________________
+void AliJetESDFillUnitArrayEMCalDigits::GetCalibrationParameters()
+{
+  // Set calibration parameters:
+  // if calibration database exists, they are read from database,
+  // otherwise, they are taken from digitizer.
+  //
+  // It is a user responsilibity to open CDB before reconstruction,
+  // for example:
+  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
+
+  //Check if calibration is stored in data base
+
+  if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
+    {
+      AliCDBEntry *entry = (AliCDBEntry*)
+        AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+      if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
+    }
+
+  if(!fCalibData)
+    printf("************* Calibration parameters not found in CDB! ****************");
+//    AliFatal("Calibration parameters not found in CDB!");
+
+
+}
+
+//____________________________________________________________________________
+Float_t  AliJetESDFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
+{
+
+  // Convert digitized amplitude into energy.
+  // Calibration parameters are taken from calibration data base for raw data,
+  // or from digitizer parameters for simulated data.
+
+  if(fCalibData){
+
+    if (fGeom==0)
+      printf("************* Did not get geometry from EMCALLoader ***************");
+//      AliFatal("Did not get geometry from EMCALLoader") ;
+
+    Int_t iSupMod = -1;
+    Int_t nModule  = -1;
+    Int_t nIphi   = -1;
+    Int_t nIeta   = -1;
+    Int_t iphi    = -1;
+    Int_t ieta    = -1;
+
+    Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
+    if(!bCell) {
+      // fGeom->PrintGeometry();
+      Error("Calibrate()"," Wrong cell id number : %i", AbsId);
+      assert(0);
+    }
+
+    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
+
+    fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
+    fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
+
+   return -fADCpedestalECA + amp * fADCchannelECA ;
+
+  }
+  else //Return energy with default parameters if calibration is not available
+    return -fADCpedestalECA + amp * fADCchannelECA ;
+
+}
+*/
+
diff --git a/JETAN/AliJetESDFillUnitArrayEMCalDigits.h b/JETAN/AliJetESDFillUnitArrayEMCalDigits.h
new file mode 100755 (executable)
index 0000000..230adba
--- /dev/null
@@ -0,0 +1,77 @@
+#ifndef ALIJETESDFILLUNITARRAYEMCALDIGITS_H
+#define ALIJETESDFILLUNITARRAYEMCALDIGITS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Jet Fill Unit Array 
+// Called by ESD Reader for jet analysis
+// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr)
+//---------------------------------------------------------------------
+
+#include "AliJetFillUnitArray.h"
+
+class AliESDCaloCluster;
+class AliJetReader;
+class AliJetESDReader;
+class AliJetESDReaderHeader;
+//class AliEMCALCalibData ;
+
+class AliJetESDFillUnitArrayEMCalDigits : public AliJetFillUnitArray
+{
+ public: 
+  AliJetESDFillUnitArrayEMCalDigits();
+  AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *fESD);
+  virtual ~AliJetESDFillUnitArrayEMCalDigits();
+  
+  // Setter
+  void SetApplyElectronCorrection(Int_t flag = 1)     {fApplyElectronCorrection = flag;}
+  void SetApplyFractionHadronicCorrection(Bool_t val) {fApplyFractionHadronicCorrection = val;}
+  void SetFractionHadronicCorrection(Double_t val)    {fFractionHadronicCorrection = val;}
+  void SetESD(AliESDEvent* const esd)                 {fESD = esd;}
+  void SetInitMult(Int_t mult)                        {fNDigitEmcal = mult;}
+  void SetInitMultCut(Int_t multcut)                  {fNDigitEmcalCut = multcut;}
+
+  // Getter
+  Int_t         GetMult()    const {return fNDigitEmcal;}
+  Int_t         GetMultCut() const {return fNDigitEmcalCut;}
+
+  // Other
+  void Exec(Option_t* const option);
+  // For calibration 
+  //  virtual Float_t Calibrate(Int_t amp, Int_t cellId) ;  // Tranforms Amp to energy
+
+ protected:
+  AliESDEvent  *fESD;                 // ESD
+  Int_t        fNIn;                  // Number of Array filled in UnitArray
+  Int_t        fOpt;                  // Detector to be used for jet reconstruction
+  Int_t        fCluster;              // Use all cells or cells in clusters for jet finding 
+  Int_t        fDebug;                // Debug option
+  Int_t        fNCEMCAL;              // Number of clusters in EMCAL
+  Int_t        fNCPHOS;               // Number of clusters in PHOS
+  Int_t        fNCCalo;               // Number of cluster in EMCAL + PHOS calorimeters
+
+  Bool_t       fApplyElectronCorrection;          // Electron correction flag
+  Bool_t       fApplyFractionHadronicCorrection;  // Fraction hadronic correction flag
+  Bool_t       fFractionHadronicCorrection;       // Fraction hadronic correction 
+
+  AliESDCaloCluster *fClus;           //! 
+  Int_t              fNDigitEmcal;    //!
+  Int_t              fNDigitEmcalCut; //!
+  // Calibration parameters... to be replaced by database
+/*   AliEMCALCalibData *fCalibData;      //! Calibration database if aval */
+/*   Float_t            fADCchannelECA;  // width of one ADC channel for EC section (GeV) */
+/*   Float_t            fADCpedestalECA; // pedestal of ADC for EC section (GeV) */
+
+
+ private:
+  AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det);
+  AliJetESDFillUnitArrayEMCalDigits &operator=(const AliJetESDFillUnitArrayEMCalDigits &det);
+
+//  void    GetCalibrationParameters(void) ;
+  
+  ClassDef(AliJetESDFillUnitArrayEMCalDigits,1) // Fill Unit Array with tpc and/or emcal information
+};
+
+#endif
diff --git a/JETAN/AliJetESDFillUnitArrayTracks.cxx b/JETAN/AliJetESDFillUnitArrayTracks.cxx
new file mode 100644 (file)
index 0000000..837fc33
--- /dev/null
@@ -0,0 +1,624 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+
+//======================================================================
+// ***July 2006
+// Fill Unit Array class 
+// Class used by AliJetESDReader to fill a UnitArray from the information extracted 
+// from the particle tracks
+// Author: magali.estienne@ires.in2p3.fr
+//======================================================================
+
+
+// --- ROOT system ---
+#include <TVector3.h>
+#include <TProcessID.h>
+
+// --- AliRoot header files ---
+#include "AliJetUnitArray.h"
+#include "AliJetESDFillUnitArrayTracks.h"
+
+// --- ROOT system ---
+class TSystem;
+class TLorentzVector;
+class TGeoManager;
+
+// --- AliRoot header files ---
+class AliJetFinder;
+class AliJetReader;
+class AliJetESDReader;
+
+ClassImp(AliJetESDFillUnitArrayTracks)
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks()
+  : AliJetFillUnitArray(),
+    fNumUnits(0),
+    fHadCorr(0),
+    fApplyMIPCorrection(kTRUE),
+    fESD(0x0),
+    fGrid0(0x0),
+    fGrid1(0x0),
+    fGrid2(0x0),
+    fGrid3(0x0),
+    fGrid4(0x0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(AliESDEvent* esd)
+  : AliJetFillUnitArray(),
+    fNumUnits(0),
+    fHadCorr(0),
+    fApplyMIPCorrection(kTRUE),
+    fESD(esd),
+    fGrid0(0x0),
+    fGrid1(0x0),
+    fGrid2(0x0),
+    fGrid3(0x0),
+    fGrid4(0x0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(const AliJetESDFillUnitArrayTracks &det)
+  : AliJetFillUnitArray(det),
+    fNumUnits(det.fNumUnits),
+    fHadCorr(det.fHadCorr),
+    fApplyMIPCorrection(det.fApplyMIPCorrection),
+    fESD(det.fESD),
+    fGrid0(det.fGrid0),
+    fGrid1(det.fGrid1),
+    fGrid2(det.fGrid2),
+    fGrid3(det.fGrid3),
+    fGrid4(det.fGrid4)
+{
+  // Copy constructor
+}
+
+//------------------------------------------------------------------
+AliJetESDFillUnitArrayTracks& AliJetESDFillUnitArrayTracks::operator=(const AliJetESDFillUnitArrayTracks& other)
+{
+  // Assignment
+
+    fNumUnits = other.fNumUnits;
+    fHadCorr = other.fHadCorr;
+    fApplyMIPCorrection = other.fApplyMIPCorrection;
+    fESD = other.fESD;
+    fGrid0 = other.fGrid0;
+    fGrid1 = other.fGrid1;
+    fGrid2 = other.fGrid2;
+    fGrid3 = other.fGrid3;
+    fGrid4 = other.fGrid4;
+
+    return (*this);
+}
+
+//____________________________________________________________________________
+void AliJetESDFillUnitArrayTracks::InitParameters()
+{
+  fHadCorr        = 0;                 // For hadron correction
+  fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
+  cout << "In AliJetESDFillUnitArrayTracks:InitParameters(), Ncells : " << fNumUnits << endl;
+
+  fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, 
+                       fPhiMax,fEtaMin,fEtaMax);
+  fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc, 
+                       fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
+
+  fIndex = fTPCGrid->GetIndexObject();
+
+  if(fDebug>20){
+    for(Int_t i=0; i<fNphi+1; i++)
+      for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
+         (*fIndex)(i,j) << endl; }
+  } 
+  if(fDebug>1) printf("\n Parameters initiated ! \n");
+}
+
+//_____________________________________________________________________________
+AliJetESDFillUnitArrayTracks::~AliJetESDFillUnitArrayTracks()
+{
+  // destructor
+}
+
+//_____________________________________________________________________________
+void AliJetESDFillUnitArrayTracks::Exec(Option_t* const /*option*/)
+{
+  //
+  // Main method.
+  // Fill the unit array with the charged particle information in ESD
+  //
+
+  fDebug = fReaderHeader->GetDebug();
+  
+  // Set parameters
+  InitParameters();
+
+  // get number of tracks in event (for the loop)
+  Int_t goodTrack = 0;
+  Int_t nt = 0;
+//(not used ?)  Int_t nt2 = 0;
+  Int_t nmax = 0;
+  Float_t pt,eta,phi;
+  //  Int_t sflag = 0;
+  TVector3 p3;
+
+  if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl;
+
+  nt = fESD->GetNumberOfTracks();
+  if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
+  // 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();  
+  fOpt = fReaderHeader->GetDetector();
+  fDZ  = fReaderHeader->GetDZ();
+
+  Int_t nTracksEmcal      = 0;
+  Int_t nTracksEmcalDZ    = 0;
+  Int_t nTracksTpc        = 0;
+  Int_t nTracksTpcOnly    = 0;
+  Int_t nTracksEmcalCut   = 0;
+  Int_t nTracksEmcalDZCut = 0;
+  Int_t nTracksTpcCut     = 0;
+  Int_t nTracksTpcOnlyCut = 0;
+
+//(not used ?)  Int_t nElements  = fTPCGrid->GetNEntries();
+//(not used ?)  Int_t nElements2 = fEMCalGrid->GetNEntries();
+  fGrid = fTPCGrid->GetGridType();
+
+//(not used ?)  Int_t nEntries = fUnitArray->GetEntries();
+
+//(not used ?)  Int_t nRefEnt = fRefArray->GetEntries();
+
+  //loop over tracks
+    nmax = nt;  
+    for (Int_t it = 0; it < nmax; it++) {
+      AliESDtrack *track;
+      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();
+
+      Float_t mass = 0.;
+      mass = track->GetMass();
+
+      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();
+       phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
+   
+      if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
+
+      sflag[goodTrack]=0;
+      if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+      cflag[goodTrack]=0;
+      if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
+
+      if(fGrid==0)
+       {
+         // Only TPC filled from its grid in its total acceptance
+         
+         Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+         Bool_t ok = kFALSE;
+         Bool_t ref = kFALSE;
+
+         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
+         TRefArray *reference = uArray->GetUnitTrackRef();
+         reference->Add(track);
+
+         Float_t unitEnergy = 0.;
+         unitEnergy = uArray->GetUnitEnergy();
+         // nTracksTpcOnly is necessary to count the number of candidate cells
+         // but it doesn't give the real multiplicity -> it will be extracted 
+         // in the jet finder through the reference to tracks
+         if(unitEnergy==0.){
+           nTracksTpcOnly++;
+           ok = kTRUE;
+            ref = kTRUE;
+         }
+
+         // Fill energy in TPC acceptance
+         uArray->SetUnitEnergy(unitEnergy + pt);
+         uArray->SetUnitMass(mass);
+
+         // Pt cut flag
+         if(uArray->GetUnitEnergy()<ptMin){
+           uArray->SetUnitCutFlag(kPtSmaller);
+         }
+         else {
+           uArray->SetUnitCutFlag(kPtHigher);
+           if(ok) nTracksTpcOnlyCut++;
+         }
+
+         // Signal flag
+         if(sflag[goodTrack] == 1) { 
+           uArray->SetUnitSignalFlag(kGood);
+            uArray->SetUnitSignalFlagC(kFALSE,kGood);
+          } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
+
+         if(uArray->GetUnitEnergy()>0 && ref){
+           if(!fProcId) {
+             fProcId = kTRUE;
+             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+           }
+           fRefArray->Add(uArray);
+         }
+       }
+    
+      if(fGrid==1)
+       {
+         Int_t nElements = fTPCGrid->GetNEntries();
+         // Fill track information in EMCAL acceptance
+           if((eta >= fEtaMin && eta <= fEtaMax) &&
+            (phi >= fPhiMin && phi <= fPhiMax))// &&
+           {
+
+             // Include dead-zones
+             if(fDZ)
+               {
+                 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);
+                 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
+                 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
+                 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
+                 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
+                 Int_t n0 = fGrid0->GetNEntries();
+                 Int_t n1 = fGrid1->GetNEntries();
+                 Int_t n2 = fGrid2->GetNEntries();
+                 Int_t n3 = fGrid3->GetNEntries();
+//(not used ?)           Int_t n4 = fGrid4->GetNEntries();
+                 
+                 if(phi >= phimin0 && phi <= phimax0){
+                   Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
+                   AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
+                   TRefArray *reference0 = uArray0->GetUnitTrackRef();
+                   reference0->Add(track);
+
+                   Float_t uEnergy0 = uArray0->GetUnitEnergy();
+                   Bool_t ok0 = kFALSE;
+                   Bool_t ref0 = kFALSE;
+                   if(uEnergy0==0.){
+                     nTracksEmcalDZ++;
+                     ok0 = kTRUE;
+                     ref0 = kTRUE;
+                   }
+                   uArray0->SetUnitEnergy(uEnergy0+pt);
+
+                   if(uArray0->GetUnitEnergy()<ptMin)
+                     uArray0->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray0->SetUnitCutFlag(kPtHigher);
+                     if(ok0) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray0->SetUnitSignalFlag(kGood);
+                      uArray0->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray0->GetUnitEnergy()>0 && ref0){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray0);
+                   }
+                 }
+
+                 if(phi >= phimin1 && phi <= phimax1){
+                   Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
+                   AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
+                   TRefArray *reference1 = uArray1->GetUnitTrackRef();
+                   reference1->Add(track);
+
+                   Float_t uEnergy1 = uArray1->GetUnitEnergy();
+                   Bool_t ok1 = kFALSE;
+                   Bool_t ref1 = kFALSE;
+                   if(uEnergy1==0.){
+                     nTracksEmcalDZ++;
+                     ok1 = kTRUE;
+                     ref1 = kTRUE;
+                   }
+                   uArray1->SetUnitEnergy(uEnergy1+pt);
+
+                   if(uArray1->GetUnitEnergy()<ptMin)
+                     uArray1->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray1->SetUnitCutFlag(kPtHigher);
+                     if(ok1) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray1->SetUnitSignalFlag(kGood);
+                      uArray1->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray1->GetUnitEnergy()>0 && ref1){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray1);
+                   }
+                 }
+
+                 if(phi >= phimin2 && phi <= phimax2){
+                   Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
+                   AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
+                   TRefArray *reference2 = uArray2->GetUnitTrackRef();
+                   reference2->Add(track);
+
+                   Float_t uEnergy2 = uArray2->GetUnitEnergy();
+                   Bool_t ok2 = kFALSE;
+                   Bool_t ref2 = kFALSE;
+                   if(uEnergy2==0.){
+                     nTracksEmcalDZ++;
+                     ok2 = kTRUE;
+                     ref2 = kTRUE;
+                   }
+                   uArray2->SetUnitEnergy(uEnergy2+pt);
+
+                   if(uArray2->GetUnitEnergy()<ptMin)
+                     uArray2->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray2->SetUnitCutFlag(kPtHigher);
+                     if(ok2) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray2->SetUnitSignalFlag(kGood);
+                      uArray2->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray2->GetUnitEnergy()>0 && ref2){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray2);
+                   }
+                 }
+
+                 if(phi >= phimin3 && phi <= phimax3){
+                   Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
+                   AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
+                   TRefArray *reference3 = uArray3->GetUnitTrackRef();
+                   reference3->Add(track);
+
+                   Float_t uEnergy3 = uArray3->GetUnitEnergy();
+                   Bool_t ok3 = kFALSE;
+                   Bool_t ref3 = kFALSE;
+                   if(uEnergy3==0.){
+                     nTracksEmcalDZ++;
+                     ok3 = kTRUE;
+                     ref3 = kTRUE;
+                   }
+                   uArray3->SetUnitEnergy(uEnergy3+pt);
+
+                   if(uArray3->GetUnitEnergy()<ptMin)
+                     uArray3->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray3->SetUnitCutFlag(kPtHigher);
+                     if(ok3) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray3->SetUnitSignalFlag(kGood);
+                      uArray3->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray3->GetUnitEnergy()>0 && ref3){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray3);
+                   }
+                 }
+
+                 if(phi >= phimin4 && phi <= phimax4){
+                   Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
+                   AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
+                   TRefArray *reference4 = uArray4->GetUnitTrackRef();
+                   reference4->Add(track);
+
+                   Float_t uEnergy4 = uArray4->GetUnitEnergy();
+                   Bool_t ok4 = kFALSE;
+                   Bool_t ref4 = kFALSE;
+                   if(uEnergy4==0.){
+                     nTracksEmcalDZ++;
+                     ok4 = kTRUE;
+                     ref4 = kTRUE;
+                   }
+                   uArray4->SetUnitEnergy(uEnergy4+pt);
+
+                   if(uArray4->GetUnitEnergy()<ptMin)
+                     uArray4->SetUnitCutFlag(kPtSmaller);
+                   else {
+                     uArray4->SetUnitCutFlag(kPtHigher);
+                     if(ok4) nTracksEmcalDZCut++;
+                   }
+                   if(sflag[goodTrack] == 1) {
+                     uArray4->SetUnitSignalFlag(kGood);
+                      uArray4->SetUnitSignalFlagC(kFALSE,kGood);
+                    } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
+
+                   if(uArray4->GetUnitEnergy()>0 && ref4){
+                     if(!fProcId){
+                       new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
+                       fProcId = kTRUE;
+                     }
+                     fRefArray->Add(uArray4);
+                   }
+                 }
+               } // end fDZ
+           
+             Int_t towerID = 0;
+             fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
+             if(towerID==-1) continue;
+             
+             AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
+             TRefArray *reference = uArray->GetUnitTrackRef();
+             reference->Add(track);
+
+             Float_t unitEnergy = uArray->GetUnitEnergy(); 
+             Bool_t ok = kFALSE;
+             Bool_t ref = kFALSE;
+             if(unitEnergy==0.){
+               nTracksEmcal++;
+               ok=kTRUE;
+               ref=kTRUE;
+             }
+
+             // Do Hadron Correction
+             // This is under construction !!!!!!!!!!!!!!!!!!!!!!!
+             // Parametrization to be added
+             if (fApplyMIPCorrection != 0) 
+               { 
+//               Float_t   hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
+//               unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
+               } //end Hadron Correction loop
+
+             uArray->SetUnitEnergy(unitEnergy + pt);
+
+             // Put a pt cut flag
+             if(uArray->GetUnitEnergy()<ptMin){
+               uArray->SetUnitCutFlag(kPtSmaller);
+             }
+             else {
+               uArray->SetUnitCutFlag(kPtHigher);
+               if(ok) nTracksEmcalCut++;
+             }
+
+             // Signal flag
+             if(sflag[goodTrack] == 1) { 
+               uArray->SetUnitSignalFlag(kGood);
+                uArray->SetUnitSignalFlagC(kFALSE,kGood);
+              } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
+
+           
+             if(uArray->GetUnitEnergy()>0 && ref){
+               if(!fProcId){
+                 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+                 fProcId = kTRUE;
+               }
+               fRefArray->Add(uArray);
+             }
+           } // end loop on EMCal acceptance cut + tracks quality
+         else{ 
+           // Outside EMCal acceptance
+         
+           Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+
+           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
+           TRefArray *reference = uArray->GetUnitTrackRef();
+           reference->Add(track);
+
+           Float_t unitEnergy2 = uArray->GetUnitEnergy(); 
+           Bool_t okout = kFALSE;
+           Bool_t refout = kFALSE;
+           if(unitEnergy2==0.){
+             nTracksTpc++;
+             okout=kTRUE;
+             refout=kTRUE;
+           }
+           // Fill energy outside emcal acceptance
+           uArray->SetUnitEnergy(unitEnergy2 + pt);
+         
+           // Pt cut flag
+           if(uArray->GetUnitEnergy()<ptMin){
+             uArray->SetUnitCutFlag(kPtSmaller);
+           }
+           else {
+             uArray->SetUnitCutFlag(kPtHigher);
+             if(okout) nTracksTpcCut++;
+           }
+           // Signal flag
+           if(sflag[goodTrack] == 1) {
+             uArray->SetUnitSignalFlag(kGood);
+              uArray->SetUnitSignalFlagC(kFALSE,kGood);
+            } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
+
+
+           if(uArray->GetUnitEnergy()>0 && refout){
+             if(!fProcId){
+               new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
+               fProcId = kTRUE;
+             }
+             fRefArray->Add(uArray);
+           }
+         }
+       } // end fGrid==1
+
+      goodTrack++;
+
+    } // end loop on entries (tpc tracks)
+
+  if(fDebug>0)
+    {
+      cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
+      cout << "goodTracks: " << goodTrack << endl;
+    }
+
+  delete[] sflag;
+  delete[] cflag;
+      
+  if(fGrid==0) {
+    fNTracks = nTracksTpcOnly;
+    fNTracksCut = nTracksTpcOnlyCut;
+    if(fDebug>10){
+      cout << "fNTracks : " << fNTracks << endl;
+      cout << "fNTracksCut : " << fNTracksCut << endl;
+    }
+  }
+  if(fGrid==1) {
+    fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
+    fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
+    if(fDebug>10){
+      cout << "fNTracks : " << fNTracks << endl;
+      cout << "fNTracksCut : " << fNTracksCut << endl;
+    }
+  }  
+  
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/JETAN/AliJetESDFillUnitArrayTracks.h b/JETAN/AliJetESDFillUnitArrayTracks.h
new file mode 100644 (file)
index 0000000..200773b
--- /dev/null
@@ -0,0 +1,64 @@
+#ifndef ALIJETESDFILLUNITARRAYTRACKS_H
+#define ALIJETESDFILLUNITARRAYTRACKS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Jet Fill Unit Array 
+// Called by ESD Reader for jet analysis
+// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr)
+//---------------------------------------------------------------------
+
+#include "AliJetFillUnitArray.h"
+#include "AliJetESDReaderHeader.h"
+
+class AliEMCALGeometry;
+class AliJetReader;
+class AliJetESDReader;
+
+class AliJetESDFillUnitArrayTracks : public AliJetFillUnitArray
+{
+ public: 
+  AliJetESDFillUnitArrayTracks();
+  AliJetESDFillUnitArrayTracks(AliESDEvent *fESD);
+  virtual ~AliJetESDFillUnitArrayTracks();
+  
+  // Setter
+  void SetHadCorrector(AliJetHadronCorrection* const corr) {fHadCorr = corr;}
+  void SetApplyMIPCorrection(Bool_t const val) {fApplyMIPCorrection = val;}
+  void SetESD(AliESDEvent* const esd) {fESD = esd;}
+  void SetGrid0(AliJetGrid* const grid0)                   {fGrid0 = grid0;}
+  void SetGrid1(AliJetGrid* const grid1)                   {fGrid1 = grid1;}
+  void SetGrid2(AliJetGrid* const grid2)                   {fGrid2 = grid2;}
+  void SetGrid3(AliJetGrid* const grid3)                   {fGrid3 = grid3;}
+  void SetGrid4(AliJetGrid* const grid4)                   {fGrid4 = grid4;}
+
+  // Getter
+  Int_t         GetHadCorrection() const {return fApplyMIPCorrection;}
+  Int_t         GetMult() const {return fNTracks;}
+  Int_t         GetMultCut() const {return fNTracksCut;}
+
+  // Other
+  void Exec(Option_t* const option);
+
+ protected:
+  Int_t                       fNumUnits;           // Number of units in the unit object array (same as num towers in EMCAL)
+  AliJetHadronCorrection      *fHadCorr;           // Pointer to Hadron Correction Object
+  Bool_t                      fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection
+  AliESDEvent                 *fESD;               // ESD
+  AliJetGrid                  *fGrid0;             // Grid used for dead zones definition
+  AliJetGrid                  *fGrid1;             // Grid used for dead zones definition
+  AliJetGrid                  *fGrid2;             // Grid used for dead zones definition
+  AliJetGrid                  *fGrid3;             // Grid used for dead zones definition
+  AliJetGrid                  *fGrid4;             // Grid used for dead zones definition
+
+ private:
+  AliJetESDFillUnitArrayTracks(const AliJetESDFillUnitArrayTracks &det);
+  AliJetESDFillUnitArrayTracks &operator=(const AliJetESDFillUnitArrayTracks &det);
+  void InitParameters();
+
+  ClassDef(AliJetESDFillUnitArrayTracks,1) 
+};
+
+#endif
index 6675a5d..6db1134 100755 (executable)
@@ -45,8 +45,8 @@
 #include "AliESD.h"
 #include "AliESDtrack.h"
 #include "AliJetDummyGeo.h"
-#include "AliJetFillUnitArrayTracks.h"
-#include "AliJetFillUnitArrayEMCalDigits.h"
+#include "AliJetESDFillUnitArrayTracks.h"
+#include "AliJetESDFillUnitArrayEMCalDigits.h"
 #include "AliJetUnitArray.h"
 #include "AliAnalysisTask.h"
 
@@ -67,9 +67,11 @@ AliJetESDReader::AliJetESDReader():
   fGrid3(0),
   fGrid4(0),
   fPtCut(0),
-  fHCorrection(0),
-  fECorrection(0),
+  fApplyElectronCorrection(kFALSE),
   fEFlag(kFALSE),
+  fApplyMIPCorrection(kTRUE),
+  fApplyFractionHadronicCorrection(kFALSE),
+  fFractionHadronicCorrection(0.3),
   fNumUnits(0),
   fDebug(0),
   fMass(0),
@@ -151,6 +153,68 @@ void AliJetESDReader::OpenInputFiles()
 
 }
 
+//__________________________________________________________
+void AliJetESDReader::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 AliJetESDReader::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 AliJetESDReader::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 AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) {
     // Connect the tree
@@ -247,7 +311,7 @@ void AliJetESDReader::CreateTasks(TChain* tree)
   // 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 = new AliJetESDFillUnitArrayTracks(); 
   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
   fFillUAFromTracks->SetGeom(fGeom);
   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
@@ -260,15 +324,17 @@ void AliJetESDReader::CreateTasks(TChain* tree)
       fFillUAFromTracks->SetGrid3(fGrid3);
       fFillUAFromTracks->SetGrid4(fGrid4);
     }
-  fFillUAFromTracks->SetHadCorrection(fHCorrection);
+  fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
   fFillUAFromTracks->SetHadCorrector(fHadCorr);
   // Create a task for to fill the neutral particle information 
-  fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits();
+  fFillUAFromEMCalDigits = new AliJetESDFillUnitArrayEMCalDigits();
   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
   fFillUAFromEMCalDigits->SetGeom(fGeom);
   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
-  fFillUAFromEMCalDigits->SetEleCorrection(fECorrection);
+  fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
+  fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
+  fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
   // Add the task to global task
   fFillUnitArray->Add(fFillUAFromTracks);
   fFillUnitArray->Add(fFillUAFromEMCalDigits);
@@ -289,7 +355,6 @@ Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
 
   fProcId = procid;
   fRefArray = refArray;
-  vector<Float_t> vtmp(3);
 
   // clear momentum array
   ClearArray();
@@ -446,15 +511,6 @@ void AliJetESDReader::InitUnitArray()
          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;
 
        }
index d42595c..e4d6b54 100755 (executable)
@@ -14,8 +14,6 @@
 // Author : magali.estienne@subatech.in2p3.fr
 //---------------------------------------------------------------------
 
-#include <vector>
-
 #include "AliJetReader.h"
 #include "AliJetUnitArray.h"
 #include "AliJetGrid.h"
@@ -50,16 +48,18 @@ class AliJetESDReader : public AliJetReader
   void       SetTPCGrid(AliJetGrid *grid)   {fTpcGrid = grid;}
   void       SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;}
   // Correction of hadronic energy
-  void       SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;}
-  void       SetHadronCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;}
-  void       SetElectronCorrection(Int_t flag = 1) {fECorrection = flag; fEFlag=kTRUE;}
+  void       SetHadronCorrector(AliJetHadronCorrection* corr) {fHadCorr = corr;}
+  void       SetApplyElectronCorrection(Int_t flag = 1) {fECorrection = flag; fEFlag=kTRUE;}
+  void       SetApplyMIPCorrection(Bool_t val);
+  void       SetApplyFractionHadronicCorrection(Bool_t val);
+  void       SetFractionHadronicCorrection(Double_t val);
 
  protected:
   AliJetDummyGeo             *fGeom;             //! EMCAL Geometry 
   TChain                     *fChain;            //! chain for reconstructed tracks
-  TChain                     *fTree;             //! tree for reconstructed tracks
+  TChain                     *fTree;             //! tree to get the EMCal geometry
   AliESDEvent                *fESD;              //! pointer to esd
-  AliJetHadronCorrectionv1   *fHadCorr;          //! Pointer to Hadron Correction Object 
+  AliJetHadronCorrection     *fHadCorr;          //! Pointer to Hadron Correction Object 
   AliJetGrid                 *fTpcGrid;          //! Pointer to grid object
   AliJetGrid                 *fEmcalGrid;        //! Pointer to grid object
   AliJetGrid                 *fGrid0;            // Pointer to grid object
@@ -68,9 +68,11 @@ class AliJetESDReader : public AliJetReader
   AliJetGrid                 *fGrid3;            // Pointer to grid object
   AliJetGrid                 *fGrid4;            // Pointer to grid object
   Float_t                     fPtCut;            // Pt cut for tracks to minimise background contribution
-  Int_t                       fHCorrection;      // Hadron correction flag
-  Int_t                       fECorrection;      // Electron correction flag
+  Int_t                       fApplyElectronCorrection;      // Electron correction flag
   Bool_t                      fEFlag;            // if (fEFlag == kFALSE) => fECorrection automatically setted
+  Bool_t                      fApplyMIPCorrection; // Apply MIP or not ? Exclusive with fApplyFractionHadronicCorrection
+  Bool_t                      fApplyFractionHadronicCorrection; // Another type of charged particle energy deposition in EMC
+  Double_t                    fFractionHadronicCorrection; // Fraction of momentum of the TPC track to be subtracted from EMC tower
   Int_t                       fNumUnits;         // Number of units in the unit object array
                                                  // (same as num towers in EMCAL)
   Int_t                       fDebug;            //! Debug option
diff --git a/JETAN/AliJetFillUnitArray.cxx b/JETAN/AliJetFillUnitArray.cxx
new file mode 100644 (file)
index 0000000..e2a91ab
--- /dev/null
@@ -0,0 +1,143 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+
+//======================================================================
+// ***July 2009
+// Fill Unit Array class 
+// Base class used by AliJetESDReader to fill a UnitArray from the information extracted 
+// from the particle tracks or emcal cells
+// Author: magali.estienne@subatech.in2p3.fr
+//======================================================================
+
+
+#include "AliJetFillUnitArray.h"
+
+// --- ROOT system ---
+class TSystem;
+class TLorentzVector;
+class TVector3;
+class TGeoManager;
+class TProcessID;
+
+// --- AliRoot header files ---
+class AliJetFinder;
+class AliJetReader;
+class AliJetESDReader;
+class AliJetESDReaderHeader;
+class AliJetUnitArray;
+
+ClassImp(AliJetFillUnitArray)
+
+//_____________________________________________________________________________
+AliJetFillUnitArray::AliJetFillUnitArray()
+  : TTask("AliJetFillUnitArray","Fill Unit Array with tpc/its and emcal information"),
+    fNTracks(0),
+    fNTracksCut(0),
+    fOpt(0),
+    fDZ(0),
+    fDebug(0),
+    fReaderHeader(0x0),
+    fMomentumArray(0x0),
+    fUnitArray(0x0),
+    fRefArray(0x0),
+    fProcId(kFALSE),
+    fTPCGrid(0x0),
+    fEMCalGrid(0x0),
+    fGeom(0x0),
+    fNphi(0),
+    fNeta(0),
+    fGrid(0),
+    fPhi2(0),
+    fEta2(0),
+    fIndex(0x0),
+    fParams(0x0),
+    fPhiMin(0),
+    fPhiMax(0),
+    fEtaMin(0),
+    fEtaMax(0),
+    fEtaBinInTPCAcc(0),
+    fPhiBinInTPCAcc(0),
+    fEtaBinInEMCalAcc(0),
+    fPhiBinInEMCalAcc(0),
+    fNbinPhi(0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetFillUnitArray::~AliJetFillUnitArray()
+{
+  // destructor
+}
+
+//_____________________________________________________________________________
+void AliJetFillUnitArray::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
+{
+  // Get the eta,phi position from the index
+
+  for(Int_t j=0; j<fNphi+1; j++) {
+    for(Int_t i=0; i<fNeta+1; i++) {
+
+      // TPC grid only 
+      //-------------------------------------
+      if(fGrid==0) {   
+       if(j*(fNeta+1)+i == index) {
+         eta = fEta2->At(i); 
+         phi = fPhi2->At(j);
+       }
+      }
+
+      // TPC-EMCAL grid
+      //-------------------------------------
+      Int_t ii = 0;
+      if(i==0) ii = 0;
+      if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i; 
+      if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
+
+      if(fGrid==1) {
+       if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
+         eta = fEta2->At(i);
+         phi = fPhi2->At(j);
+       }  
+
+       if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) && 
+          ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
+         if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
+         else eta = fEta2->At(i);
+         phi = fPhi2->At(j);
+       }
+
+       if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
+         eta = fEta2->At(i);
+         phi = fPhi2->At(j);
+       }
+      }
+    }
+  }
+}
+
+//_____________________________________________________________________________
+Float_t  AliJetFillUnitArray::EtaToTheta(Float_t arg)
+{
+  return 2.*atan(exp(-arg));
+}
+
+
+
+
+
+
+
diff --git a/JETAN/AliJetFillUnitArray.h b/JETAN/AliJetFillUnitArray.h
new file mode 100644 (file)
index 0000000..1c190b5
--- /dev/null
@@ -0,0 +1,122 @@
+#ifndef ALIJETFILLUNITARRAY_H
+#define ALIJETFILLUNITARRAY_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// Base class used to fill Unit Array 
+// Called by ESD Reader for jet analysis
+// Author: Magali Estienne (magali.estienne@subatech.in2p3.fr)
+//---------------------------------------------------------------------
+
+//#include <Riostream.h>
+#include <TMath.h>
+#include <TTask.h>
+
+#include <TArrayD.h>
+#include <TRefArray.h>
+#include <TClonesArray.h>
+
+#include "AliJetReaderHeader.h"
+#include "AliJetAODReaderHeader.h"
+#include "AliJetDummyGeo.h"
+#include "AliJetGrid.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliJetHadronCorrection.h"
+
+class Riostream;
+
+class AliEMCALGeometry;
+class AliJetReader;
+class AliJetESDReader;
+
+class AliJetFillUnitArray : public TTask
+{
+ public: 
+  AliJetFillUnitArray();
+  virtual ~AliJetFillUnitArray();
+  AliJetFillUnitArray(const AliJetFillUnitArray &det);
+  AliJetFillUnitArray &operator=(const AliJetFillUnitArray &det);
+  
+  // Setter
+  virtual void SetReaderHeader(AliJetReaderHeader* const readerHeader) {fReaderHeader = readerHeader;}
+  virtual void SetGeom(AliJetDummyGeo* const geom)                     {fGeom = geom;}
+  virtual void SetMomentumArray(TClonesArray* const momentumArray)     {fMomentumArray = momentumArray;}
+  virtual void SetUnitArray(TClonesArray* const unitArray)             {fUnitArray = unitArray;}
+  virtual void SetRefArray(TRefArray* const refArray)                  {fRefArray = refArray;}
+  virtual void SetTPCGrid(AliJetGrid* const grid)                      {fTPCGrid = grid;}
+  virtual void SetEMCalGrid(AliJetGrid* const grid)                    {fEMCalGrid = grid;}
+  virtual void SetProcId(Bool_t id)                              {fProcId = id;}
+  virtual void SetGrid0(AliJetGrid */*grid0*/)                    {;}
+  virtual void SetGrid1(AliJetGrid */*grid1*/)                    {;}
+  virtual void SetGrid2(AliJetGrid */*grid2*/)                    {;}
+  virtual void SetGrid3(AliJetGrid */*grid3*/)                    {;}
+  virtual void SetGrid4(AliJetGrid */*grid4*/)                    {;}
+  virtual void SetHadCorrector(AliJetHadronCorrection* /*corr*/)  {;}
+  virtual void SetApplyMIPCorrection(Bool_t /*val*/)              {;}
+  virtual void SetESD(AliESDEvent */*esd*/)                       {;}
+  virtual void SetAOD(AliAODEvent */*aod*/)                       {;}
+  virtual void SetApplyElectronCorrection(Int_t /*flag*/)         {;}
+  virtual void SetApplyFractionHadronicCorrection(Bool_t /*val*/) {;}
+  virtual void SetFractionHadronicCorrection(Double_t /*val*/)    {;}
+  virtual void SetInitMult(Int_t /*mult*/)                        {;}
+  virtual void SetInitMultCut(Int_t /*multcut*/)                  {;}
+
+  // Getter
+  virtual TClonesArray* GetUnitArray() const {return fUnitArray;}
+  virtual TRefArray*    GetRefArray() const {return fRefArray;}
+  virtual void          GetEtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi);
+  virtual Int_t         GetNeta() const {return fNeta;}
+  virtual Int_t         GetNphi() const {return fNphi;}
+  virtual Int_t         GetMult() const {return 0;} // To be checked
+  virtual Int_t         GetMultCut() const {return 0;} // To be checked
+  virtual Bool_t        GetProcId() const {return fProcId;}
+
+  // Other
+  virtual void          Exec(Option_t* const /*option*/) {;}
+  virtual Float_t       EtaToTheta(Float_t arg);
+  virtual void          InitParameters() {;}
+
+ protected:
+  Int_t                 fNTracks;         // Number of tracks stored in UnitArray
+  Int_t                 fNTracksCut;      // Number of tracks stored in UnitArray with a pt cut 
+  Int_t                 fOpt;             // Detector to be used for jet reconstruction
+  Bool_t                fDZ;              // Use or not dead zones
+  Int_t                 fDebug;           // Debug option
+
+  AliJetReaderHeader   *fReaderHeader;    // ReaderHeader
+  TClonesArray         *fMomentumArray;   // MomentumArray
+  TClonesArray         *fUnitArray;       // UnitArray
+  TRefArray            *fRefArray;        // UnitArray
+  Bool_t                fProcId;          // Bool_t for TProcessID synchronization  
+  AliJetGrid           *fTPCGrid;         // Define filled grid
+  AliJetGrid           *fEMCalGrid;       // Define filled grid
+  AliJetDummyGeo       *fGeom;            // Define EMCal geometry
+
+  Int_t                 fNphi;            // number of points in the grid:   phi
+  Int_t                 fNeta;            //               "                 eta
+  Int_t                 fGrid;            // Select the grid acceptance you want to fill
+                                          // 0 = TPC acceptance, 1 = TPC-EMCal acceptance
+  TArrayD              *fPhi2;            // grid points in phi
+  TArrayD              *fEta2;            // grid points in eta
+  TMatrixD             *fIndex;           // grid points in (phi,eta) 
+  TMatrixD             *fParams;          // matrix of parameters in the grid points  
+  Float_t               fPhiMin;          // EMCal acceptance
+  Float_t               fPhiMax;          // EMCal acceptance
+  Float_t               fEtaMin;          // EMCal acceptance
+  Float_t               fEtaMax;          // EMCal acceptance
+  Int_t                 fEtaBinInTPCAcc;  // Number of bins in Eta in TPC acceptance
+  Int_t                 fPhiBinInTPCAcc;  // Number of bins in phi in TPC acceptance
+  Int_t                 fEtaBinInEMCalAcc;// Number of bins in Eta in EMCal acceptance
+  Int_t                 fPhiBinInEMCalAcc;// Number of bins in phi in EMCal acceptance
+  Int_t                 fNbinPhi;         // Number of phi bins
+
+ private:
+
+
+  ClassDef(AliJetFillUnitArray,1) // Fill Unit Array with tpc and/or emcal information
+};
+
+#endif
index 2523989..bbdaa67 100644 (file)
 
 #include <Riostream.h>
 #include <TFile.h>
-#include <TClonesArray.h>
-#include <TProcessID.h>
 
 #include "AliJetFinder.h"
 #include "AliJet.h"
 #include "AliAODJet.h"
-#include "AliJetReader.h"
-#include "AliJetReaderHeader.h"
 #include "AliJetControlPlots.h"
 #include "AliLeading.h"
 #include "AliAODEvent.h"
 #include "AliJetUnitArray.h"
 
+class TProcessID;
+class TClonesArray;
+
 ClassImp(AliJetFinder)
 
 AliJetFinder::AliJetFinder():
@@ -252,7 +251,6 @@ Bool_t AliJetFinder::ProcessEvent2()
   FindJets();
   
   Int_t         nEntRef    = ref->GetEntries();
-  vector<Float_t> vtmp;
 
   for(Int_t i=0; i<nEntRef; i++)
     { 
@@ -261,10 +259,11 @@ Bool_t AliJetFinder::ProcessEvent2()
       ((AliJetUnitArray*)ref->At(i))->SetUnitEnergy(0.);
       ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag(kPtSmaller);
       ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag2(kPtSmaller);
-      ((AliJetUnitArray*)ref->At(i))->SetUnitPxPyPz(kTRUE,vtmp);
       ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlag(kBad);
+      ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlagC(kTRUE,kBad);
       ((AliJetUnitArray*)ref->At(i))->SetUnitDetectorFlag(kTpc);
       ((AliJetUnitArray*)ref->At(i))->SetUnitFlag(kOutJet);
+      ((AliJetUnitArray*)ref->At(i))->ClearUnitTrackRef();
 
       // Reset process ID
       AliJetUnitArray* uA = (AliJetUnitArray*)ref->At(i);
index a782d7c..3d6a100 100755 (executable)
 //          magali.estienne@subatech.in2p3.fr
 //---------------------------------------------------------------------
 
+//#include <Riostream.h>
+//#include <vector>
+
 #include <TObject.h>
 #include "AliAODJet.h"
 #include "AliJetReader.h"
+#include "AliJetHeader.h"
+#include "AliJet.h"
+#include "AliJetReaderHeader.h"
 
 class TFile;
 class TChain;
-class AliJet;
-class AliJetHeader;
+//class AliJet;
+//class AliJetHeader;
 class AliJetControlPlots;
 class AliLeading;
 class AliAODJet;
@@ -53,7 +59,7 @@ class AliJetFinder : public TObject
   virtual void          Reset() {fNAODjets = 0;}
   virtual void          FindJets() {}
   virtual void          FindJetsC(){}
-  virtual void          WriteJHeaderToFile() { }
+  virtual void          WriteJHeaderToFile() {}
   // some methods to allow steering from the outside
   virtual Bool_t        ProcessEvent();
   virtual Bool_t        ProcessEvent2();
@@ -63,7 +69,7 @@ class AliJetFinder : public TObject
   virtual void          ConnectAODNonStd(AliAODEvent* aod,const char* bname);
   virtual TTree         *MakeTreeJ(char* name);
   virtual void          WriteHeaders();
-  virtual void          WriteJetsToFile() {;}
+  virtual void          WriteJetsToFile() const {;}
   virtual void          TestJet() {;}
 
  protected:
index 2364b97..88a3a4c 100755 (executable)
@@ -29,8 +29,8 @@
 #include "AliJetReaderHeader.h"
 #include "AliESDEvent.h"
 #include "AliHeader.h"
-#include "AliJetFillUnitArrayTracks.h" 
-#include "AliJetFillUnitArrayEMCalDigits.h"
+#include "AliJetESDFillUnitArrayTracks.h" 
+#include "AliJetESDFillUnitArrayEMCalDigits.h"
 #include "AliJetUnitArray.h"
 #include "AliJetHadronCorrectionv1.h"
 
@@ -51,10 +51,9 @@ AliJetReader::AliJetReader():
   fSignalFlag(0),
   fCutFlag(0),
   fUnitArray(new TClonesArray("AliJetUnitArray",60000)),
-  fUnitArrayNoCuts(new TClonesArray("AliJetUnitArray",60000)),
   fArrayInitialised(0),
-  fFillUAFromTracks(new AliJetFillUnitArrayTracks()), 
-  fFillUAFromEMCalDigits(new AliJetFillUnitArrayEMCalDigits()),
+  fFillUAFromTracks(new AliJetESDFillUnitArrayTracks()), 
+  fFillUAFromEMCalDigits(new AliJetESDFillUnitArrayEMCalDigits()),
   fNumCandidate(0),
   fNumCandidateCut(0),
   fHadronCorrector(0),
@@ -82,11 +81,6 @@ AliJetReader::~AliJetReader()
       delete fUnitArray;
   }
   
-  if (fUnitArrayNoCuts) {
-    fUnitArrayNoCuts->Delete();
-    delete fUnitArrayNoCuts;
-  }
-
   if (fFillUnitArray) {
     delete fFillUnitArray;
   }
index 457f388..72aa1d4 100755 (executable)
@@ -22,10 +22,9 @@ class AliJetReaderHeader;
 class AliESDEvent;
 class AliHeader;
 class AliJetUnitArray;
-class AliJetHadronCorrectionv1;
+class AliJetHadronCorrection;
 class AliJet;
-class AliJetFillUnitArrayTracks;
-class AliJetFillUnitArrayEMCalDigits;
+class AliJetFillUnitArray;
 
 class AliJetReader : public TObject 
 {
@@ -37,8 +36,6 @@ class AliJetReader : public TObject
   virtual TClonesArray*       GetMomentumArray()     const {return fMomentumArray;}
   virtual TRefArray*          GetReferences()        const {return 0;}   
   virtual TClonesArray        *GetUnitArray() const {return fUnitArray;}  
-  virtual TClonesArray        *GetUnitArrayNoCuts()  const {return fUnitArrayNoCuts;} 
-
   virtual AliJetReaderHeader* GetReaderHeader()      const {return fReaderHeader;}
   virtual AliHeader           *GetAliHeader() const  {return fAliHeader;}
   virtual Int_t               GetSignalFlag(Int_t i) const {return fSignalFlag[i];}
@@ -50,17 +47,16 @@ class AliJetReader : public TObject
   // Setters
   virtual Bool_t FillMomentumArray() {return kTRUE;}
   virtual Bool_t ReadEventLoader(Int_t) {return kTRUE;}
-  virtual void   FillUnitArrayFromTPCTracks(Int_t) {}     // temporarily not used
-  virtual void   FillUnitArrayFromEMCALHits() {}          // temporarily not used
-  virtual void   FillUnitArrayFromEMCALDigits(Int_t) {}   // temporarily not used
   virtual void   InitUnitArray() {}
   virtual void   InitParameters() {}
   virtual void   CreateTasks(TChain* /*tree*/) {}
   virtual Bool_t ExecTasks(Bool_t /*procid*/, TRefArray* /*refArray*/) {return kFALSE;}
   // Correction of hadronic energy 
-  virtual void   SetHadronCorrector(AliJetHadronCorrectionv1*) {;} 
-  virtual void   SetHadronCorrection(Int_t) {;} 
-  virtual void   SetElectronCorrection(Int_t) {;} 
+  virtual void   SetHadronCorrector(AliJetHadronCorrection*) {;} 
+  virtual void   SetApplyMIPCorrection(Bool_t /*val*/){;}
+  virtual void   SetApplyFractionHadronicCorrection(Bool_t /*val*/){;}
+  virtual void   SetFractionHadronicCorrection(Double_t /*val*/){;}
+  virtual void   SetApplyElectronCorrection(Int_t /*flag*/) {;}
   virtual void   SetReaderHeader(AliJetReaderHeader* header) 
   {fReaderHeader = header;}
   virtual void   SetESD(AliESDEvent* esd) { fESD = esd;}
@@ -89,14 +85,13 @@ class AliJetReader : public TObject
                                                            // from the underlying event
   TArrayI                          fCutFlag;               // to flag if a particle passed the pt cut or not
   TClonesArray                    *fUnitArray;             // array of digit position and energy 
-  TClonesArray                    *fUnitArrayNoCuts;       //|| array of digit position and energy 
   Bool_t                           fArrayInitialised;      // To check that array of units is initialised  
-  AliJetFillUnitArrayTracks       *fFillUAFromTracks;      // For charged particle task
-  AliJetFillUnitArrayEMCalDigits  *fFillUAFromEMCalDigits; // For neutral particle task
+  AliJetFillUnitArray             *fFillUAFromTracks;      // For charged particle task
+  AliJetFillUnitArray             *fFillUAFromEMCalDigits; // For neutral particle task
   Int_t                            fNumCandidate;          // Number of entries different from zero in unitarray
   Int_t                            fNumCandidateCut;       // Number of entries different from zero in unitarray
                                                            // which pass pt cut
-  AliJetHadronCorrectionv1        *fHadronCorrector;       //! Pointer to hadronic correction 
+  AliJetHadronCorrection          *fHadronCorrector;       //! Pointer to hadronic correction 
   Int_t                            fHCorrection;           //  Hadron correction flag 
   Int_t                            fECorrection;           //  Electron correction flag 
   Bool_t                           fEFlag;                 //  Electron correction flag 
index 5b55049..28ffcd6 100644 (file)
 //*-- Author: Sarah Blyth (LBL/UCT)
 // --
 // Revised Version for JETAN 
-// -- Magali Estienne (IReS)
-
-//#include <vector>
-
-#include <TVector3.h>
-#include <TLorentzVector.h>
-#include <TClonesArray.h>
+// -- Magali Estienne (magali.estienne@subatech.in2p3.fr)
 
 #include "AliJetUnitArray.h"
 
+class TVector3;
+class TLorentzVector;
+class TClonesArray;
+
 ClassImp(AliJetUnitArray)
 
 AliJetUnitArray::AliJetUnitArray():
@@ -53,10 +51,11 @@ AliJetUnitArray::AliJetUnitArray():
   fUnitPy(0.),
   fUnitPz(0.),
   fUnitMass(0.),
-  fV(0),
   fVc(0),
   fVn(0),
-  fVet(0)
+  fUnitTrackRef(new TRefArray),
+  fUnitCellRef(new TRefArray),
+  fUnitClusterRef(new TRefArray)
 {
   // Default constructor
 }  
@@ -74,16 +73,17 @@ AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t
   fUnitFlag(inout),
   fUnitCutFlag(cut),
   fUnitCutFlag2(cut2),
-  fUnitSignalFlag(signal), 
+  fUnitSignalFlag(kBad), 
   fUnitDetectorFlag(det),
   fUnitPx(0.),
   fUnitPy(0.),
   fUnitPz(0.),
   fUnitMass(mass),
-  fV(0),
   fVc(0),
   fVn(0),
-  fVet(0)
+  fUnitTrackRef(new TRefArray),
+  fUnitCellRef(new TRefArray),
+  fUnitClusterRef(new TRefArray)
 {
   //abs ID (in a eta,phi grid, track ID in ESD, eta, phi, energy, px, py, pz, Deta, Dphi, detector flag, in/out jet, mass
 
@@ -109,18 +109,94 @@ AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t
   fUnitPy(py),
   fUnitPz(pz),
   fUnitMass(mass),
-  fV(0),
   fVc(0),
   fVn(0),
-  fVet(0)
+  fUnitTrackRef(new TRefArray),
+  fUnitCellRef(new TRefArray),
+  fUnitClusterRef(new TRefArray)
 {
   // Constructor 2
 }
  
+AliJetUnitArray::AliJetUnitArray(const AliJetUnitArray& rUnit):
+  TObject(rUnit),
+  fUnitEnergy(rUnit.fUnitEnergy),
+  fUnitEta(rUnit.fUnitEta),
+  fUnitPhi(rUnit.fUnitPhi),
+  fUnitDeta(rUnit.fUnitDeta),
+  fUnitDphi(rUnit.fUnitDphi),
+  fUnitID(rUnit.fUnitID),
+  fUnitTrackID(rUnit.fUnitTrackID),
+  fUnitNum(rUnit.fUnitNum),
+  fUnitClusterID(rUnit.fUnitClusterID),
+  fUnitFlag(rUnit.fUnitFlag),
+  fUnitCutFlag(rUnit.fUnitCutFlag),
+  fUnitCutFlag2(rUnit.fUnitCutFlag2),
+  fUnitSignalFlag(rUnit.fUnitSignalFlag),
+  fUnitDetectorFlag(rUnit.fUnitDetectorFlag),
+  fUnitPx(rUnit.fUnitPx),
+  fUnitPy(rUnit.fUnitPy),
+  fUnitPz(rUnit.fUnitPz),
+  fUnitMass(rUnit.fUnitMass),
+  fVc(rUnit.fVc),
+  fVn(rUnit.fVn),
+  fUnitTrackRef(rUnit.fUnitTrackRef),
+  fUnitCellRef(rUnit.fUnitCellRef),
+  fUnitClusterRef(rUnit.fUnitClusterRef)
+{
+  // Copy constructor
+}
+
+AliJetUnitArray& AliJetUnitArray::operator=(const AliJetUnitArray& rhs)
+{
+  // Assignment
+  fUnitEnergy = rhs.fUnitEnergy;
+  fUnitEta = rhs.fUnitEta;
+  fUnitPhi = rhs.fUnitPhi;
+  fUnitDeta = rhs.fUnitDeta;
+  fUnitDphi = rhs.fUnitDphi;
+  fUnitID = rhs.fUnitID;
+  fUnitTrackID = rhs.fUnitTrackID;
+  fUnitNum = rhs.fUnitNum;
+  fUnitClusterID = rhs.fUnitClusterID;
+  fUnitFlag = rhs.fUnitFlag;
+  fUnitCutFlag = rhs.fUnitCutFlag;
+  fUnitCutFlag2 = rhs.fUnitCutFlag2;
+  fUnitSignalFlag = rhs.fUnitSignalFlag;
+  fUnitDetectorFlag = rhs.fUnitDetectorFlag;
+  fUnitPx = rhs.fUnitPx;
+  fUnitPy = rhs.fUnitPy;
+  fUnitPz = rhs.fUnitPz;
+  fUnitMass = rhs.fUnitMass;
+  fVc = rhs.fVc;
+  fVn = rhs.fVn;
+  fUnitTrackRef = rhs.fUnitTrackRef;
+  fUnitCellRef = rhs.fUnitCellRef;
+  fUnitClusterRef = rhs.fUnitClusterRef;
+
+  return *this;
+
+}
+
+
 //------------------------------------------------------------------------
 AliJetUnitArray::~AliJetUnitArray()
 {
   // Destructor 
+  delete fUnitTrackRef;
+  delete fUnitCellRef;
+  delete fUnitClusterRef;
+
+}
+
+void AliJetUnitArray::ClearUnitTrackRef()
+{
+  fUnitTrackRef->Clear();
+}
+
+void AliJetUnitArray::ClearUnitCellRef()
+{
+  fUnitCellRef->Clear();
 }
 
 //------------------------------------------------------------------------
@@ -145,32 +221,6 @@ void AliJetUnitArray::SetUnitSignalFlagN(Bool_t init, AliJetFinderUnitSignalFlag
   else fVn.push_back(flag);
 }
 
-//------------------------------------------------------------------------
-void AliJetUnitArray::SetUnitEtN(Bool_t init, Float_t et)
-{
-  // Set transverse energy of the neutral cell
-  if(init){
-    if(!fVet.empty())
-      fVet.clear();
-  }
-  else fVet.push_back(et);
-}
-
-
-//------------------------------------------------------------------------
-void AliJetUnitArray::SetUnitPxPyPz(Bool_t init, vector<Float_t> v3)
-{
-  // Set momentum components of the charged particle
-  if(init)
-    {
-      if(!fV.empty()){
-       fV.clear();
-      }
-    }
-  else{
-    fV.push_back(v3);
-  }
-}
 
 //------------------------------------------------------------------------
 Bool_t AliJetUnitArray::GetUnitSignalFlagC(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagc)
@@ -197,47 +247,6 @@ Bool_t AliJetUnitArray::GetUnitSignalFlagN(Int_t ind, AliJetFinderUnitSignalFlag
 }
 
 //------------------------------------------------------------------------
-Bool_t AliJetUnitArray::GetUnitEtN(Int_t ind, Float_t &et)
-{
-  // Get transverse energy of the neutral cell
-  if(ind <= (Int_t)fVet.size())
-    {
-      et = (Float_t)fVet[ind];
-      return kTRUE;
-    }
-  else return kFALSE;
-}
-
-//------------------------------------------------------------------------
-Bool_t AliJetUnitArray::GetUnitPxPyPz(Int_t ind, Float_t &px, Float_t &py, Float_t &pz)
-{
-  // Get momentum components of the charged particle
-  if(ind <= (Int_t)fV.size())
-    {
-      px = (Float_t)fV[ind][0];
-      py = (Float_t)fV[ind][1];
-      pz = (Float_t)fV[ind][2];
-      return kTRUE;
-    }
-  else return kFALSE;
-}
-
-//------------------------------------------------------------------------
-Bool_t AliJetUnitArray::GetUnitPxPyPzE(Int_t ind, Float_t &px, Float_t &py, Float_t &pz, Float_t &en)
-{
-// Get 4-momentum components of the charged particle
-  if(ind <= (Int_t)fV.size())
-    {
-      px = (Float_t)fV[ind][0];
-      py = (Float_t)fV[ind][1];
-      pz = (Float_t)fV[ind][2];
-      en = TMath::Sqrt(px*px+py*py+pz*pz);
-      return kTRUE;
-    }
-  else return kFALSE;
-}
-
-//------------------------------------------------------------------------
 Float_t  AliJetUnitArray::EtaToTheta(Float_t arg) const
 {
   // Eta to theta transformation
index 1673cda..5300bc2 100644 (file)
 //
 //
 
-#include <Riostream.h>
-#include <vector>
-
+#include <Riostream.h> 
 #include <TObject.h>
-#include <TList.h>
-#include <TVector3.h>
+#include <TRefArray.h>
+#include <TMath.h>
+#include "AliESDtrack.h"
+#include "AliAODTrack.h"
+#include "AliVParticle.h"
+#include "AliESDCaloCells.h"
+#include "AliAODCaloCells.h"
+#include "AliESDCaloCluster.h"
+#include "AliAODCaloCluster.h"
 #include "AliJetFinderTypes.h"
 
+class TList;
+class TVector3;
+
 class AliJetUnitArray : public TObject
 {
  public:
@@ -68,36 +76,36 @@ class AliJetUnitArray : public TObject
   { 
     fUnitDetectorFlag = detectorflag; 
   } 
-  void SetUnitPxPyPz(Bool_t init, vector<Float_t> v3);
-  void SetUnitEtN(Bool_t init, Float_t et); // Added for background studies
   void SetUnitSignalFlagC(Bool_t init, AliJetFinderUnitSignalFlagType_t flag);
   void SetUnitSignalFlagN(Bool_t init, AliJetFinderUnitSignalFlagType_t flag);
   void SetUnitMass(Float_t mass) {fUnitMass = mass;}
+  void SetUnitTrackRef(TRefArray* trackref)     {fUnitTrackRef = trackref;}
+  void SetUnitCellRef(TRefArray* cellref)       {fUnitCellRef = cellref;}
+  void SetUnitClusterRef(TRefArray* clusterref) {fUnitClusterRef = clusterref;}
 
   // Getter
-  Float_t GetUnitEnergy() const            {return fUnitEnergy;}
-  Float_t GetUnitEta() const               {return fUnitEta;}
-  Float_t GetUnitPhi() const               {return fUnitPhi;} 
-  Float_t GetUnitPx() const                {return fUnitPx;}
-  Float_t GetUnitPy() const                {return fUnitPy;}
-  Float_t GetUnitPz() const                {return fUnitPz;}
-  Float_t GetUnitDeta() const              {return fUnitDeta;}
-  Float_t GetUnitDphi() const              {return fUnitDphi;}         
-  Int_t   GetUnitID() const                {return fUnitID;}
-  Int_t   GetUnitTrackID() const           {return fUnitTrackID;}
-  Int_t   GetUnitEntries() const           {return fUnitNum;}
-  Int_t   GetUnitClusterID() const         {return fUnitClusterID;}
-  Float_t GetUnitMass() const              {return fUnitMass;}
-  Bool_t  GetUnitPxPyPz(Int_t ind, Float_t &px, Float_t &py, Float_t &pz);
-  Bool_t  GetUnitPxPyPzE(Int_t ind, Float_t &px, Float_t &py, Float_t &pz, Float_t &en);
-  Bool_t  GetUnitEtN(Int_t ind, Float_t &et); // Added for background studies
-  Bool_t  GetUnitSignalFlagC(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagc);
-  Bool_t  GetUnitSignalFlagN(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagn);
-  Int_t   GetUnitVectorSize() {return fV.size();}
-  Int_t   GetUnitVectorSizeC() {return fVc.size();}
-  Int_t   GetUnitVectorSizeN() {return fVn.size();}
-
-  Float_t EtaToTheta(Float_t arg) const;
+  Float_t     GetUnitEnergy() const            {return fUnitEnergy;}
+  Float_t     GetUnitEta() const               {return fUnitEta;}
+  Float_t     GetUnitPhi() const               {return fUnitPhi;} 
+  Float_t     GetUnitPx() const                {return fUnitPx;}
+  Float_t     GetUnitPy() const                {return fUnitPy;}
+  Float_t     GetUnitPz() const                {return fUnitPz;}
+  Float_t     GetUnitDeta() const              {return fUnitDeta;}
+  Float_t     GetUnitDphi() const              {return fUnitDphi;}         
+  Int_t       GetUnitID() const                {return fUnitID;}
+  Int_t       GetUnitTrackID() const           {return fUnitTrackID;}
+  Int_t       GetUnitEntries() const           {return fUnitNum;}
+  Int_t       GetUnitClusterID() const         {return fUnitClusterID;}
+  Float_t     GetUnitMass() const              {return fUnitMass;}
+  Bool_t      GetUnitSignalFlagC(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagc);
+  Bool_t      GetUnitSignalFlagN(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagn);
+  TRefArray*  GetUnitTrackRef() const          {return fUnitTrackRef;}
+  TRefArray*  GetUnitCellRef() const           {return fUnitCellRef;}
+  TRefArray*  GetUnitClusterRef() const        {return fUnitClusterRef;}
+
+  Float_t     EtaToTheta(Float_t arg) const;
+  void        ClearUnitTrackRef();
+  void        ClearUnitCellRef();
 
   AliJetFinderUnitFlagType_t GetUnitFlag() const     
   {
@@ -121,6 +129,8 @@ class AliJetUnitArray : public TObject
   }
 
  protected:
+  AliJetUnitArray(const AliJetUnitArray& rUnit);
+  AliJetUnitArray& operator = (const AliJetUnitArray& rhs);
   Bool_t operator>  ( AliJetUnitArray unit1) const;
   Bool_t operator<  ( AliJetUnitArray unit1) const;
   Bool_t operator== ( AliJetUnitArray unit1) const;
@@ -143,10 +153,11 @@ class AliJetUnitArray : public TObject
   Float_t                            fUnitPy;           // Py of charged track
   Float_t                            fUnitPz;           // Pz of charged track
   Float_t                            fUnitMass;         // Mass of particle
-  vector< vector< Float_t > >                  fV;      //|| vector to store part information in each cell
   vector< AliJetFinderUnitSignalFlagType_t >   fVc;     //|| added for background studies
   vector< AliJetFinderUnitSignalFlagType_t >   fVn;     //|| added for background studies
-  vector< Float_t >                            fVet;    //|| added for background studies          
+  TRefArray*                         fUnitTrackRef;     //! pointer to array of references to esd tracks
+  TRefArray*                         fUnitCellRef;      //! pointer to array of references to esd cells
+  TRefArray*                         fUnitClusterRef;   //! pointer to array of references to esd clusters
 
   ClassDef(AliJetUnitArray,1)
 
index 9e11e17..04decb9 100644 (file)
 // Author: magali.estienne@subatech.in2p3.fr
 //---------------------------------------------------------------------
 
-#include <Riostream.h>
-#include <vector>
-
-#include <TArrayF.h>
 #include <TClonesArray.h>
-#include <TFile.h>
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TLorentzVector.h>
 #include <TMath.h>
 #include <TRefArray.h>
+#include "TFile.h"
 
 #include "AliUA1JetFinderV2.h"
 #include "AliUA1JetHeaderV1.h"
 #include "AliJetUnitArray.h"
-#include "AliJetReaderHeader.h"
-#include "AliJetReader.h"
-#include "AliJet.h"
-#include "AliAODJet.h"
 
+class TArrayF;
+class TFile;
+class AliJetReader;
+class AliAODJet;
 
 ClassImp(AliUA1JetFinderV2)
 
@@ -109,8 +105,8 @@ void AliUA1JetFinderV2::FindJetsC()
     ptT[i]  = lv->Pt();
     etaT[i] = lv->Eta();
     phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
-    cFlagT[i] = fReader->GetCutFlag(i); // Temporarily added
-    sFlagT[i] = fReader->GetSignalFlag(i); // Temporarily added peut-etre a mettre apres cut en pt !!!
+    cFlagT[i] = fReader->GetCutFlag(i); 
+    sFlagT[i] = fReader->GetSignalFlag(i); 
     
     if (fReader->GetCutFlag(i) != 1) continue;
     fLego->Fill(etaT[i], phiT[i], ptT[i]);
@@ -127,10 +123,10 @@ void AliUA1JetFinderV2::FindJetsC()
   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
   
   // arrays to hold jets
-  Float_t* etaJet    = new Float_t[30];
-  Float_t* phiJet    = new Float_t[30];
-  Float_t* etJet     = new Float_t[30];
-  Float_t* etsigJet  = new Float_t[30]; //signal et in jet
+  Float_t* etaJet    = new Float_t[30];  // eta jet
+  Float_t* phiJet    = new Float_t[30];  // phi jet
+  Float_t* etJet     = new Float_t[30];  // et jet
+  Float_t* etsigJet  = new Float_t[30];  // signal et in jet
   Float_t* etallJet  = new Float_t[30];  // total et in jet (tmp variable)
   Int_t*   ncellsJet = new Int_t[30];
   Int_t*   multJet   = new Int_t[30];
@@ -138,7 +134,7 @@ void AliUA1JetFinderV2::FindJetsC()
   Float_t* etaJetOk    = new Float_t[30];
   Float_t* phiJetOk    = new Float_t[30];
   Float_t* etJetOk     = new Float_t[30];
-  Float_t* etsigJetOk  = new Float_t[30]; //signal et in jet
+  Float_t* etsigJetOk  = new Float_t[30];  // signal et in jet
   Float_t* etallJetOk  = new Float_t[30];  // total et in jet (tmp variable)
   Int_t*   ncellsJetOk = new Int_t[30];
   Int_t*   multJetOk   = new Int_t[30];
@@ -181,7 +177,6 @@ void AliUA1JetFinderV2::FindJetsC()
     //subtract background
     Float_t etbgTotalN = 0.0; //new background
     if(header->GetBackgMode() == 1) // standard
-      //      SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
       SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
     if(header->GetBackgMode() == 2) //cone
       SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
@@ -210,6 +205,7 @@ void AliUA1JetFinderV2::FindJetsC()
     phiJetOk[p]    = phiJet[idx[p]];
     etJetOk[p]     = etJet[idx[p]];
     etallJetOk[p]  = etJet[idx[p]];
+    etsigJetOk[p]  = etsigJet[idx[p]];
     ncellsJetOk[p] = ncellsJet[idx[p]];
     multJetOk[p]   = multJet[idx[p]];
   }
@@ -325,8 +321,9 @@ void AliUA1JetFinderV2::FindJets()
   Int_t              nCand    = fReader->GetNumCandidate();
   Int_t              nCandCut = fReader->GetNumCandidateCut();
   Int_t              nIn      = fUnit->GetEntries();
-  Float_t            fPtMin   = fReader->GetReaderHeader()->GetPtCut();
+  Float_t            ptMin   = fReader->GetReaderHeader()->GetPtCut();
 
+  fDebug   = fReader->GetReaderHeader()->GetDebug();
   if (nIn == 0) return;
 
   Int_t nCandidateCut = 0;
@@ -337,36 +334,34 @@ void AliUA1JetFinderV2::FindJets()
 
   // local arrays for input No Cuts
   // Both pt < ptMin and pt > ptMin
-  Float_t* ptT       = new Float_t[nCandidate];
-  Float_t* en2T      = new Float_t[nCandidate];
-  Float_t* pt2T      = new Float_t[nCandidate];
-  Int_t*   detT      = new Int_t[nCandidate]; 
-  Float_t* etaT      = new Float_t[nCandidate];
-  Float_t* phiT      = new Float_t[nCandidate];
-  Float_t* cFlagT    = new Float_t[nCandidate];
-  Float_t* cFlag2T   = new Float_t[nCandidate];
-  Float_t* sFlagT    = new Float_t[nCandidate];
-  Float_t* cClusterT = new Float_t[nCandidate];
-  Int_t*   vectT     = new Int_t[nCandidate];
-  Int_t    loop1     = 0;
-  Int_t*   injet     = new Int_t[nCandidate];
-  Int_t*   sflag     = new Int_t[nCandidate];
-  vector< vector<Float_t> > pxT;
-  vector< vector<Float_t> > pyT;
-  vector< vector<Float_t> > pzT;
+  Float_t*   ptT       = new Float_t[nCandidate];
+  Float_t*   en2T      = new Float_t[nCandidate];
+  Float_t*   pt2T      = new Float_t[nCandidate];
+  Int_t*     detT      = new Int_t[nCandidate]; 
+  Float_t*   etaT      = new Float_t[nCandidate];
+  Float_t*   phiT      = new Float_t[nCandidate];
+  Float_t*   cFlagT    = new Float_t[nCandidate];
+  Float_t*   cFlag2T   = new Float_t[nCandidate];
+  Float_t*   sFlagT    = new Float_t[nCandidate];
+  Float_t*   cClusterT = new Float_t[nCandidate];
+  Int_t*     vectT     = new Int_t[nCandidate];
+  Int_t      loop1     = 0;
+  Int_t*     injet     = new Int_t[nCandidate];
+  Int_t*     sflag     = new Int_t[nCandidate];
+  TRefArray* trackRef  = new TRefArray();
 
   //total energy in array
   Float_t  etbgTotal = 0.0;
-  TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
+  TH1F*    hPtTotal  = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
 
   // Input cell info
-  Float_t *etCell = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
-  Float_t *etaCell = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
-  Float_t *phiCell = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
-  Int_t   *flagCell = new Int_t[nIn];   //! Cell phi - Extracted from UnitArray
-  Float_t *etCell2 = new Float_t[nIn];  //! Cell Energy - Extracted from UnitArray
-  Float_t *etaCell2 = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
-  Float_t *phiCell2 = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
+  Float_t *etCell    = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
+  Float_t *etaCell   = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
+  Float_t *phiCell   = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
+  Int_t   *flagCell  = new Int_t[nIn];   //! Cell phi - Extracted from UnitArray
+  Float_t *etCell2   = new Float_t[nIn];  //! Cell Energy - Extracted from UnitArray
+  Float_t *etaCell2  = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
+  Float_t *phiCell2  = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
   Int_t   *flagCell2 = new Int_t[nIn];  //! Cell phi - Extracted from UnitArray
 
   // Information extracted from fUnitArray
@@ -376,25 +371,13 @@ void AliUA1JetFinderV2::FindJets()
       // Recover particle information from UnitArray
       
       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
-      
+      TRefArray* ref = uArray->GetUnitTrackRef();
+      Int_t nRef = ref->GetEntries();
+
       if(uArray->GetUnitEnergy()>0.){
-       vector<Float_t> vtmpx; 
-       vector<Float_t> vtmpy; 
-       vector<Float_t> vtmpz; 
-       for(Int_t j=0; j<uArray->GetUnitVectorSize();j++)
-         {
-           Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
-           uArray->GetUnitPxPyPz(j,x,y,z);
-           vtmpx.push_back(x);
-           vtmpy.push_back(y);
-           vtmpz.push_back(z);
-         }
-       pxT.push_back(vtmpx);
-       pyT.push_back(vtmpy);
-       pzT.push_back(vtmpz);
-       vtmpx.clear();
-       vtmpy.clear();
-       vtmpz.clear();
+
+       for(Int_t jpart=0; jpart<nRef;jpart++)
+         trackRef->Add((AliVTrack*)ref->At(jpart));
        ptT[loop1]   = uArray->GetUnitEnergy();
         detT[loop1]  = uArray->GetUnitDetectorFlag();
        etaT[loop1]  = uArray->GetUnitEta();
@@ -402,7 +385,7 @@ void AliUA1JetFinderV2::FindJets()
        cFlagT[loop1]= uArray->GetUnitCutFlag();   // pt cut tpc
        cFlag2T[loop1]= uArray->GetUnitCutFlag2(); // pt cut emcal
        sFlagT[loop1]= uArray->GetUnitSignalFlag();
-       vectT[loop1] = uArray->GetUnitVectorSize();
+       vectT[loop1] = nRef;
        if(cFlagT[loop1] == 1 || cFlag2T[loop1] == 1) {
          pt2T[loop1] = 0.;
          en2T[loop1] = 0.;
@@ -414,11 +397,13 @@ void AliUA1JetFinderV2::FindJets()
          }
          if(detT[loop1]==0){ // TPC+ITS
            Float_t pt = 0.;
-           for(Int_t j=0; j<vectT[loop1];j++){
+           for(Int_t j=0; j<nRef;j++){
              Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
-             uArray->GetUnitPxPyPz(j,x,y,z);
+             x = ((AliVTrack*)ref->At(j))->Px();
+             y = ((AliVTrack*)ref->At(j))->Py();
+             z = ((AliVTrack*)ref->At(j))->Pz();
              pt = TMath::Sqrt(x*x+y*y);
-             if(pt>fPtMin) {
+             if(pt>ptMin) {
                pt2T[loop1] += pt;
                en2T[loop1] += pt;
                hPtTotal->Fill(pt);
@@ -430,11 +415,13 @@ void AliUA1JetFinderV2::FindJets()
            Float_t ptCTot = 0.;
            Float_t pt = 0.;
            Float_t enC = 0.;
-           for(Int_t j=0; j<vectT[loop1];j++) {
+           for(Int_t j=0; j<nRef;j++){
              Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
-             uArray->GetUnitPxPyPz(j,x,y,z);
+             x = ((AliVTrack*)ref->At(j))->Px();
+             y = ((AliVTrack*)ref->At(j))->Py();
+             z = ((AliVTrack*)ref->At(j))->Pz();
              pt = TMath::Sqrt(x*x+y*y);
-             if(pt>fPtMin) {
+             if(pt>ptMin) {
                pt2T[loop1]+=pt;
                en2T[loop1]+=pt;
                hPtTotal->Fill(pt);
@@ -466,12 +453,14 @@ void AliUA1JetFinderV2::FindJets()
         }
         if(uArray->GetUnitDetectorFlag()==0){ // TPC case
           Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
-          for(Int_t j=0; j<uArray->GetUnitVectorSize();j++)
+          for(Int_t j=0; j<nRef;j++)
             {
               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
-              uArray->GetUnitPxPyPz(j,x,y,z);
+             x = ((AliVTrack*)ref->At(j))->Px();
+             y = ((AliVTrack*)ref->At(j))->Py();
+             z = ((AliVTrack*)ref->At(j))->Pz();
               pt = TMath::Sqrt(x*x+y*y);
-              if(pt>fPtMin) {
+              if(pt>ptMin) {
                 et1 += pt;
                 et2 += pt;
               }
@@ -490,12 +479,14 @@ void AliUA1JetFinderV2::FindJets()
           Float_t ptCTot = 0.;
           Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
           Float_t enC = 0.;
-          for(Int_t j=0; j<uArray->GetUnitVectorSize();j++)
+          for(Int_t j=0; j<nRef;j++)
             {
               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
-              uArray->GetUnitPxPyPz(j,x,y,z);
+             x = ((AliVTrack*)ref->At(j))->Px();
+             y = ((AliVTrack*)ref->At(j))->Py();
+             z = ((AliVTrack*)ref->At(j))->Pz();
               pt = TMath::Sqrt(x*x+y*y);
-              if(pt>fPtMin) {
+              if(pt>ptMin) {
                 et1 += pt;
                 et2 += pt;
               }
@@ -625,6 +616,7 @@ void AliUA1JetFinderV2::FindJets()
       phiJetOk[p]    = phiJet[idx[p]];
       etJetOk[p]     = etJet[idx[p]];
       etallJetOk[p]  = etJet[idx[p]];
+      etsigJetOk[p]  = etsigJet[idx[p]];
       ncellsJetOk[p] = ncellsJet[idx[p]];
       multJetOk[p]   = multJet[idx[p]];
     }
@@ -682,10 +674,7 @@ void AliUA1JetFinderV2::FindJets()
   fJets->SetEtaIn(etaT);
   fJets->SetPhiIn(phiT);
   fJets->SetPtIn(ptT);
-  fJets->SetVectorSizeIn(vectT);
-  fJets->SetVectorPxIn(pxT);
-  fJets->SetVectorPyIn(pyT);
-  fJets->SetVectorPzIn(pzT);
+  fJets->SetTrackReferences(trackRef);
   fJets->SetDetectorFlagIn(detT);
   fJets->SetEtAvg(etbgTotal/(2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax()-header->GetLegoPhiMin())));
 
@@ -695,9 +684,8 @@ void AliUA1JetFinderV2::FindJets()
   delete pt2T;
   delete etaT;
   delete phiT;
-  pxT.clear();
-  pyT.clear();
-  pzT.clear();
+  trackRef->Delete();
+  delete trackRef;
   delete detT;
   delete cFlagT;
   delete cFlag2T;
@@ -740,13 +728,18 @@ void AliUA1JetFinderV2::FindJets()
 
 ////////////////////////////////////////////////////////////////////////
 void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, 
-                                   Int_t* flagCell, Float_t* etCell2, Float_t* etaCell2, Float_t* phiCell2, 
-                                   Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, 
-                                   Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
-                                   Float_t* etallJet, Int_t* ncellsJet)
+                                   Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2, 
+                                   const Float_t* phiCell2, const Int_t* flagCell2, Float_t etbgTotal, 
+                                   Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet, 
+                                   Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet)
 {
-  
+  //
+  // Main method for jet finding
+  // UA1 base cone finder
+  //
+
   Int_t nCell  = nIn; 
+  fDebug = fReader->GetReaderHeader()->GetDebug();
 
   // Dump lego
   // Check enough space! *to be done*
@@ -1137,10 +1130,10 @@ void AliUA1JetFinderV2::RunAlgoritmC(Float_t etbgTotal, Double_t dEtTotal, Int_t
 }
 
 ////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, Float_t* ptT, 
-                                     Int_t*vectT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* cFlag2T, 
-                                     Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, 
-                                     Float_t* etsigJet, Int_t* multJet, Int_t* injet)
+void AliUA1JetFinderV2::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, const Float_t* ptT, 
+                                     const Int_t*vectT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, 
+                                     const Float_t* cFlag2T, const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, 
+                                     const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet)
 {
   //
   // Background subtraction using cone method but without correction in dE/deta distribution
@@ -1149,7 +1142,8 @@ void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN,
   
   //calculate energy inside and outside cones
   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
-  fOpt = fReader->GetReaderHeader()->GetDetector();
+  Int_t              fOpt = fReader->GetReaderHeader()->GetDetector();
+  fDebug = fReader->GetReaderHeader()->GetDebug();
   Float_t rc= header->GetRadius();
   Float_t etIn[30];
   Float_t etOut = 0;
@@ -1342,9 +1336,9 @@ void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN,
 }
 
 ////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::SubtractBackgC(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN,
-                                      Float_t* ptT, Float_t* etaT, Float_t* phiT,
-                                      Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
+void AliUA1JetFinderV2::SubtractBackgC(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN,
+                                      const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
+                                      Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
                                       Int_t* multJet, Int_t* injet)
 {
   //background subtraction using cone method but without correction in dE/deta distribution
@@ -1408,10 +1402,10 @@ void AliUA1JetFinderV2::SubtractBackgC(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN,
 
 
 ////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotalN,
-                     Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT,
-                      Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
-                      Int_t* multJet, Int_t* injet)
+void AliUA1JetFinderV2::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN,
+                     const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, 
+                     const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, 
+                     Float_t* etsigJet, Int_t* multJet, Int_t* injet)
 {
 
   //background subtraction using statistical method
@@ -1473,9 +1467,9 @@ void AliUA1JetFinderV2::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotal
 }
 
 ////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT,
-                                         Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT,
-                                         Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
+void AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT,
+                                         Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
+                                         Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
                                          Int_t* multJet, Int_t* injet)
 {
   // Cone background subtraction method taking into acount dEt/deta distribution
@@ -1587,9 +1581,9 @@ void AliUA1JetFinderV2::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTota
 }
 
 ////////////////////////////////////////////////////////////////////////
-void AliUA1JetFinderV2::SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN,
-                     Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT,
-                     Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
+void AliUA1JetFinderV2::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
+                     Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
+                     Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
                      Int_t* multJet, Int_t* injet)
 {
   // Ratio background subtraction method taking into acount dEt/deta distribution
index eef1429..5c64d88 100644 (file)
@@ -12,9 +12,8 @@
 // (version in c++)
 //---------------------------------------------------------------------
 
-#include <vector>
-
 #include "AliJetFinder.h"
+
 class AliUA1JetHeaderV1;
 class TH2F;
 class TChain;
@@ -34,45 +33,44 @@ class AliUA1JetFinderV2 : public AliJetFinder
                   Float_t* etallJet, Int_t* ncellsJet);
 
   void RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, 
-                  Int_t* flagCell, Float_t* etCell2, Float_t* etaCell2, Float_t* phiCell2, 
-                  Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, 
+                  Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2, const Float_t* phiCell2, 
+                  const Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, 
                   Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
                   Float_t* etallJet, Int_t* ncellsJet);
    
-  void SubtractBackgC(Int_t& nIn, Int_t&nJ, Float_t&EtbgTotalN,
-                      Float_t* ptT, Float_t* etaT, Float_t* phiT,
-                      Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
+  void SubtractBackgC(const Int_t& nIn, const Int_t&nJ, Float_t&EtbgTotalN,
+                      const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
+                      Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
                       Float_t* etsigJet,Int_t* multJet, Int_t* injet);
 
-  void SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&EtbgTotalN, Float_t* ptT, Int_t* vectT, 
-                    Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* cFlag2T, 
-                    Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, 
+  void SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&EtbgTotalN, const Float_t* ptT, const Int_t* vectT, 
+                    const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* cFlag2T, 
+                    const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, 
                     Float_t* etsigJet, Int_t* multJet, Int_t* injet);
 
-  void SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& EtbgTotalN,
-                      Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT,
-                      Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
+  void SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& EtbgTotalN,
+                      Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
+                      Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
                       Float_t* etsigJet, Int_t* multJet, Int_t* injet);
 
-  void SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& EtbgTotalN,
-                      Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT,
-                      Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
+  void SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& EtbgTotalN,
+                      Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
+                      Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
                       Float_t* etsigJet, Int_t* multJet, Int_t* injet);
 
-  void SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&EtbgTotalN,
-                      Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT,
-                      Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
+  void SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&EtbgTotalN,
+                        const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, 
+                        const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
                       Float_t* etsigJet, Int_t* multJet, Int_t* injet);
   void Reset();
   void InitTask(TChain* tree);
   void WriteJHeaderToFile();
-
  protected:
   AliUA1JetFinderV2(const AliUA1JetFinderV2& rJetF1);
   AliUA1JetFinderV2& operator = (const AliUA1JetFinderV2& rhsf);
-  TH2F           * fLego;           //Lego Histo
-  Int_t fDebug;
-  Int_t fOpt;
+  TH2F           * fLego;           // Lego Histo
+  Int_t fDebug;                     // Debug option
+  Int_t fOpt;                       // Detector option (charged only or charged+neutral)
 
   ClassDef(AliUA1JetFinderV2,1)
 };
index d870cca..5d2a01f 100644 (file)
 #pragma link C++ class AliJetHadronCorrection+; 
 #pragma link C++ class AliJetHadronCorrectionv0+;
 #pragma link C++ class AliJetHadronCorrectionv1+; 
-#pragma link C++ class AliJetFillUnitArrayTracks+;
-#pragma link C++ class AliJetFillUnitArrayEMCalDigits+;
+#pragma link C++ class AliJetFillUnitArray+;
+#pragma link C++ class AliJetESDFillUnitArrayTracks+;
+#pragma link C++ class AliJetESDFillUnitArrayEMCalDigits+;
+#pragma link C++ class AliJetAODFillUnitArrayTracks+;
+#pragma link C++ class AliJetAODFillUnitArrayEMCalDigits+;
 #pragma link C++ class AliJetDummyGeo+;
 #pragma link C++ class AliJetDummyShishKebabTrd1Module+;
 #pragma link C++ class AliAnalysisTaskJets+;
index 8ce13b1..31220ee 100644 (file)
@@ -11,7 +11,9 @@ SRCS =        AliJet.cxx AliJetHeader.cxx \
        AliCdfJetFinder.cxx AliCdfJetHeader.cxx \
        AliJetGrid.cxx AliJetUnitArray.cxx AliJetHadronCorrection.cxx \
        AliJetHadronCorrectionv0.cxx AliJetHadronCorrectionv1.cxx \
-       AliJetFillUnitArrayTracks.cxx AliJetFillUnitArrayEMCalDigits.cxx\
+        AliJetFillUnitArray.cxx        AliJetESDFillUnitArrayTracks.cxx \
+       AliJetESDFillUnitArrayEMCalDigits.cxx AliJetAODFillUnitArrayTracks.cxx \
+       AliJetAODFillUnitArrayEMCalDigits.cxx \
        AliJetDummyGeo.cxx AliJetDummyShishKebabTrd1Module.cxx\
        AliJetFinderTypes.cxx \
        AliAnalysisTaskJets.cxx \