- Use UnitArray to store track and emcal information.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 13 Sep 2007 08:12:57 +0000 (08:12 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 13 Sep 2007 08:12:57 +0000 (08:12 +0000)
- Hadron correction.
(Magali Estienne)

26 files changed:
JETAN/AliJetAODReader.cxx
JETAN/AliJetControlPlots.h
JETAN/AliJetDummyGeo.cxx
JETAN/AliJetDummyGeo.h
JETAN/AliJetESDReader.cxx
JETAN/AliJetESDReader.h
JETAN/AliJetFillUnitArrayEMCalDigits.cxx [new file with mode: 0644]
JETAN/AliJetFillUnitArrayEMCalDigits.h [new file with mode: 0644]
JETAN/AliJetFillUnitArrayTracks.cxx
JETAN/AliJetFillUnitArrayTracks.h
JETAN/AliJetFinder.cxx
JETAN/AliJetHadronCorrectionv1.cxx
JETAN/AliJetHadronCorrectionv1.h
JETAN/AliJetMCReader.cxx
JETAN/AliJetReader.cxx
JETAN/AliJetReader.h
JETAN/AliJetReaderHeader.cxx
JETAN/AliJetReaderHeader.h
JETAN/AliJetUnitArray.cxx
JETAN/AliJetUnitArray.h
JETAN/AliUA1JetFinderV2.cxx [new file with mode: 0644]
JETAN/AliUA1JetFinderV2.h [new file with mode: 0644]
JETAN/AliUA1JetHeaderV1.cxx
JETAN/AliUA1JetHeaderV1.h
JETAN/JETANLinkDef.h
JETAN/libJETAN.pkg

index 592ac61..bb28631 100644 (file)
@@ -24,6 +24,8 @@
 #include <TSystem.h>
 #include <TLorentzVector.h>
 #include <TVector3.h>
+#include <TChain.h>
+
 #include "AliJetAODReader.h"
 #include "AliJetAODReaderHeader.h"
 #include "AliAODEvent.h"
index d1e305f..48c6bda 100755 (executable)
@@ -18,6 +18,9 @@ class TFile;
 class TClonesArray;
 class TH1I;
 class TH1D;
+class TH1F;
+class TH1;
+
 class AliJetReader;
 class AliJet;
 
@@ -29,19 +32,19 @@ class AliJetControlPlots : public TObject
 
   // setter
   // getters
-  TH1I *GetNJetsH() {return fNJetsH;}
-  TH1I *GetMultH() {return fMultH;}
-  TH1D *GetPhiH() {return fPhiH;}
-  TH1D *GetFractionInJetH() {return fInJetH;}
-  TH1D *GetEneH() {return fEneH;}
-  TH1D *GetPtH() {return fPtH;}
-  TH1D *GetEtaH() {return fEtaH;}
-  TH1D *GetFragH() {return fFragH;}
-  TH1D *GetFragLnH() {return fFragLnH;}
-  TH1D *GetFragrH() {return fFragrH;}
-  TH1D *GetFragLnrH() {return fFragLnrH;}
-  TH1D *GetShapeH() {return fShapeH;}
-  TH1D *GetShaperH() {return fShaperH;}  
+  TH1I *GetNJetsH()         const {return fNJetsH;}
+  TH1I *GetMultH()          const {return fMultH;}
+  TH1D *GetPhiH()           const {return fPhiH;}
+  TH1D *GetFractionInJetH() const {return fInJetH;}
+  TH1D *GetEneH()           const {return fEneH;}
+  TH1D *GetPtH()            const {return fPtH;}
+  TH1D *GetEtaH()           const {return fEtaH;}
+  TH1D *GetFragH()          const {return fFragH;}
+  TH1D *GetFragLnH()        const {return fFragLnH;}
+  TH1D *GetFragrH()         const {return fFragrH;}
+  TH1D *GetFragLnrH()       const {return fFragLnrH;}
+  TH1D *GetShapeH()         const {return fShapeH;}
+  TH1D *GetShaperH()        const {return fShaperH;}  
   
   // others
   void FillHistos(AliJet *j);
index b2a1700..7e5010e 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
+
+//
+// Temporarily added to define part of the EMCal geometry
+// necessary for the jet finder
+//
+
+#include <assert.h>
+
+// --- Root header files ---
+#include <TVector3.h>
+#include <TGeoManager.h>
+#include <TGeoMatrix.h>
+#include <TGeoNode.h>
+
+// --- AliRoot header files ---
 #include "AliJetDummyGeo.h"
 
 ClassImp(AliJetDummyGeo)
+
+AliJetDummyGeo::AliJetDummyGeo():
+    TObject(),
+    fArm1PhiMin(80.0), 
+    fArm1PhiMax(200.0),
+    fArm1EtaMin(-0.7), 
+    fArm1EtaMax(+0.7), 
+    fNumberOfSuperModules(12), 
+    fSteelFrontThick(0.0), 
+    fIPDistance(428.0), 
+    fZLength(),
+    fPhiGapForSM(2.), 
+    fNPhi(12),       
+    fNZ(24),      
+    fPhiModuleSize(12.26 - fPhiGapForSM / Float_t(fNPhi)), 
+    fEtaModuleSize(fPhiModuleSize), 
+    fNPHIdiv(2), 
+    fNETAdiv(2), 
+    fNECLayers(77), 
+    fECScintThick(0.16), 
+    fECPbRadThickness(0.16), 
+    fSampling(12.327),
+    fTrd1Angle(1.5), 
+    fNCellsInModule(fNPHIdiv*fNETAdiv), 
+    fNCellsInSupMod(fNCellsInModule*fNPhi*fNZ), 
+    fNCells(fNCellsInSupMod*fNumberOfSuperModules-fNCellsInSupMod), 
+    fLongModuleSize(fNECLayers*(fECScintThick + fECPbRadThickness)), 
+    f2Trd1Dx2(fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.)), 
+    fShellThickness(TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2)+fSteelFrontThick),
+    fEtaMaxOfTRD1(0.67064) // Value extracted from ShishKebab
+{
+  // Constructor
+  // Local coordinates
+  fParSM[0] = GetShellThickness()/2.;        
+  fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
+  fParSM[2] = 350./2.;
+
+  fZLength    = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage
+  fEnvelop[0] = fIPDistance; // mother volume inner radius
+  fEnvelop[1] = fIPDistance + fShellThickness; // mother volume outer r.
+  fEnvelop[2] = 1.00001*fZLength; // add some padding for mother volume. 
+
+  // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006 
+  fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
+  fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
+  fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
+  fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
+  fPhiCentersOfSM[0]    = TMath::PiOver2();
+  for(int i=1; i<=4; i++) { // from 2th ro 9th
+    fPhiBoundariesOfSM[2*i]   = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
+    fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
+    fPhiCentersOfSM[i]        = fPhiCentersOfSM[0]    + 20.*TMath::DegToRad()*i;
+  }
+  fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
+  fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
+  fPhiCentersOfSM[5]     = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.; 
+
+  //  for(int i=0; i<fNumberOfSuperModules; i++) fMatrixOfSM[i] = 0;
+
+  fCentersOfCellsEtaDir.Set(fNZ*fNETAdiv);
+  fCentersOfCellsXDir.Set(fNZ*fNETAdiv);
+  fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
+  fEtaCentersOfCells.Set(fNZ*fNETAdiv*fNPhi*fNPHIdiv);
+  fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
+
+  int nphism  = GetNumberOfSuperModules()/2;
+  double dphi = (GetArm1PhiMax() - GetArm1PhiMin())/nphism;
+  double rpos = (GetEnvelop(0) + GetEnvelop(1))/2.;
+  double phi, phiRad, xpos, ypos, zpos;
+  for(int i=0; i<nphism; i++){
+    phi    = GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 90, 110, 130, 150, 170, 190
+    phiRad = phi*TMath::Pi()/180.;
+    xpos = rpos * TMath::Cos(phiRad);
+    ypos = rpos * TMath::Sin(phiRad);
+    zpos = fParSM[2];
+    if(i==5) {
+      xpos += (fParSM[1]/2. * TMath::Sin(phiRad)); 
+      ypos -= (fParSM[1]/2. * TMath::Cos(phiRad));
+    }
+    // pozitive z
+    int ind = 2*i;
+    TGeoRotation *geoRot0 = new TGeoRotation("geoRot0", 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
+    fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
+                                         xpos,ypos, zpos, geoRot0);
+    // negaive z
+    ind++;
+    double phiy = 90. + phi + 180.;
+    if(phiy>=360.) phiy -= 360.;
+    TGeoRotation *geoRot1 = new TGeoRotation("geoRot1", 90.0, phi, 90.0, phiy, 180.0, 0.0);
+    fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
+                                         xpos,ypos,-zpos, geoRot1);
+  } // for
+
+}
+
+AliJetDummyGeo::AliJetDummyGeo(const AliJetDummyGeo& geom):
+    TObject(),  
+    fArm1PhiMin(geom.fArm1PhiMin),
+    fArm1PhiMax(geom.fArm1PhiMax),
+    fArm1EtaMin(geom.fArm1EtaMin),
+    fArm1EtaMax(geom.fArm1EtaMax),
+    fNumberOfSuperModules(geom.fNumberOfSuperModules),
+    fSteelFrontThick(geom.fSteelFrontThick),
+    fIPDistance(geom.fIPDistance),
+    fPhiGapForSM(geom.fPhiGapForSM),
+    fNPhi(geom.fNPhi),
+    fNZ(geom.fNZ),
+    fPhiModuleSize(geom.fPhiModuleSize),
+    fEtaModuleSize(geom.fEtaModuleSize),
+    fNPHIdiv(geom.fNPHIdiv),
+    fNETAdiv(geom.fNETAdiv),
+    fNECLayers(geom.fNECLayers),
+    fECScintThick(geom.fECScintThick),
+    fECPbRadThickness(geom.fECPbRadThickness),
+    fSampling(geom.fSampling),
+    fTrd1Angle(geom.fTrd1Angle),
+    fNCellsInModule(geom.fNCellsInModule),
+    fNCellsInSupMod(geom.fNCellsInSupMod),
+    fNCells(geom.fNCells),
+    fLongModuleSize(geom.fLongModuleSize),
+    f2Trd1Dx2(geom.f2Trd1Dx2),
+    fShellThickness(geom.fShellThickness),
+    fEtaMaxOfTRD1(geom.fEtaMaxOfTRD1) 
+{
+  // Constructor
+  // Local coordinates
+  fParSM[0] = GetShellThickness()/2.;        
+  fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
+  fParSM[2] = 350./2.;
+
+  // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006 
+  fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
+  fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
+  fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
+  fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
+  fPhiCentersOfSM[0]    = TMath::PiOver2();
+  for(int i=1; i<=4; i++) { // from 2th ro 9th
+    fPhiBoundariesOfSM[2*i]   = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
+    fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
+    fPhiCentersOfSM[i]        = fPhiCentersOfSM[0]    + 20.*TMath::DegToRad()*i;
+  }
+  fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
+  fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
+  fPhiCentersOfSM[5]     = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.; 
+
+  //  for(int i=0; i<fNumberOfSuperModules; i++) fMatrixOfSM[i] = 0;
+
+  fCentersOfCellsEtaDir.Set(fNZ*fNETAdiv);
+  fCentersOfCellsXDir.Set(fNZ*fNETAdiv);
+  fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
+  fEtaCentersOfCells.Set(fNZ*fNETAdiv*fNPhi*fNPHIdiv);
+  fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
+
+  int nphism  = GetNumberOfSuperModules()/2;
+  double dphi = (GetArm1PhiMax() - GetArm1PhiMin())/nphism;
+  double rpos = (GetEnvelop(0) + GetEnvelop(1))/2.;
+  double phi, phiRad, xpos, ypos, zpos;
+  for(int i=0; i<nphism; i++){
+    phi    = GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 90, 110, 130, 150, 170, 190
+    phiRad = phi*TMath::Pi()/180.;
+    xpos = rpos * TMath::Cos(phiRad);
+    ypos = rpos * TMath::Sin(phiRad);
+    zpos = fParSM[2];
+    if(i==5) {
+      xpos += (fParSM[1]/2. * TMath::Sin(phiRad)); 
+      ypos -= (fParSM[1]/2. * TMath::Cos(phiRad));
+    }
+    // pozitive z
+    int ind = 2*i;
+    TGeoRotation *geoRot0 = new TGeoRotation("geoRot0", 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
+    fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
+                                         xpos,ypos, zpos, geoRot0);
+    // negaive z
+    ind++;
+    double phiy = 90. + phi + 180.;
+    if(phiy>=360.) phiy -= 360.;
+    TGeoRotation *geoRot1 = new TGeoRotation("geoRot1", 90.0, phi, 90.0, phiy, 180.0, 0.0);
+    fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
+                                         xpos,ypos,-zpos, geoRot1);
+
+    delete geoRot0;
+    delete geoRot1;
+
+  } // for
+
+}
+
+AliJetDummyGeo::~AliJetDummyGeo()
+{
+  // Destructor
+  delete [] fMatrixOfSM;
+}
+
+//------------------------------------------------------------------------------------
+void AliJetDummyGeo::EtaPhiFromIndex(Int_t absId, Float_t& eta, Float_t& phi)
+{
+  // Nov 16, 2006- float to double
+  // version for TRD1 only
+  static TVector3 vglob;
+  GetGlobal(absId, vglob);
+  eta = vglob.Eta();
+  phi = vglob.Phi();
+
+}
+
+//------------------------------------------------------------------------------------
+void AliJetDummyGeo::GetGlobal(const Double_t *loc, Double_t *glob, int ind) const
+{
+  // Figure out the global numbering of a given supermodule from the
+  // local numbering
+  // Alice numbering - Jun 03,2006
+  //  if(fMatrixOfSM[0] == 0) GetTransformationForSM();
+
+  if(ind>=0 && ind < GetNumberOfSuperModules()) {
+    fMatrixOfSM[ind]->LocalToMaster(loc, glob);
+  }
+}
+
+//------------------------------------------------------------------------------------
+void AliJetDummyGeo::GetGlobal(Int_t absId , double glob[3]) const
+{ 
+  // Alice numbering scheme - Jun 03, 2006
+  static Int_t nSupMod, nModule, nIphi, nIeta;
+  static double loc[3];
+
+  glob[0]=glob[1]=glob[2]=0.0; // bad case
+  if(RelPosCellInSModule(absId, loc)) {
+    GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
+    fMatrixOfSM[nSupMod]->LocalToMaster(loc, glob);
+  }
+}
+
+//------------------------------------------------------------------------------------
+void AliJetDummyGeo::GetGlobal(Int_t absId , TVector3 &vglob) const
+{ 
+  // Alice numbering scheme - Jun 03, 2006
+  static Double_t glob[3];
+
+  GetGlobal(absId, glob);
+  vglob.SetXYZ(glob[0], glob[1], glob[2]);
+
+}
+
+//------------------------------------------------------------------------------------
+Bool_t AliJetDummyGeo::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const
+{
+  // Alice numbering scheme - Jun 03, 2006
+  loc[0] = loc[1] = loc[2]=0.0;
+  if(RelPosCellInSModule(absId, loc[0],loc[1],loc[2])) {
+    return kTRUE;
+  }
+  return kFALSE;
+}
+
+//------------------------------------------------------------------------------------
+Bool_t AliJetDummyGeo::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const
+{
+  // Look to see what the relative position inside a given cell is
+  // for a recpoint.
+  // Alice numbering scheme - Jun 08, 2006
+  // In:
+  // absId   - cell is as in Geant,     0<= absId   < fNCells;
+  // OUT:
+  // xr,yr,zr - x,y,z coordinates of cell with absId inside SM 
+
+  // Shift index taking into account the difference between standard SM 
+  // and SM of half size in phi direction
+  const Int_t phiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
+  static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
+  if(!CheckAbsCellId(absId)) return kFALSE;
+
+  GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
+  GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta); 
+  xr = fCentersOfCellsXDir.At(ieta);
+  zr = fCentersOfCellsEtaDir.At(ieta);
+
+  if(nSupMod<10) {
+    yr = fCentersOfCellsPhiDir.At(iphi);
+  } else {
+    yr = fCentersOfCellsPhiDir.At(iphi + phiIndexShift);
+  }
+
+  return kTRUE;
+}
+
+//------------------------------------------------------------------------------------
+Bool_t  AliJetDummyGeo::CheckAbsCellId(Int_t absId) const
+{ 
+  // May 31, 2006; only trd1 now
+  if(absId<0 || absId >= fNCells) return kFALSE;
+  else                            return kTRUE;
+}
+
+//------------------------------------------------------------------------------------
+Bool_t AliJetDummyGeo::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const
+{ 
+  // 21-sep-04; 19-oct-05;
+  // May 31, 2006; ALICE numbering scheme:
+  // 
+  // In:
+  // absId   - cell is as in Geant,     0<= absId   < fNCells;
+  // Out:
+  // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
+  // nModule  - module number in SM,    0<= nModule  < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
+  // nIphi   - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv; 
+  // nIeta   - cell number in eta driection inside module; 0<= nIeta < fNETAdiv; 
+  // 
+  static Int_t tmp=0, sm10=0;
+  if(!CheckAbsCellId(absId)) return kFALSE;
+
+  sm10 = fNCellsInSupMod*10;
+  if(absId >= sm10) { // 110 degree case; last two supermodules  
+    nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
+    tmp     = (absId-sm10) % (fNCellsInSupMod/2);
+  } else {
+    nSupMod = absId / fNCellsInSupMod;
+    tmp     = absId % fNCellsInSupMod;
+  }
+
+  nModule  = tmp / fNCellsInModule;
+  tmp     = tmp % fNCellsInModule;
+  nIphi   = tmp / fNPHIdiv;
+  nIeta   = tmp % fNPHIdiv;
+
+  return kTRUE;
+}
+
+//------------------------------------------------------------------------------------
+void AliJetDummyGeo::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta, int &iphi, int &ieta) const
+{ 
+  // 
+  // Added nSupMod; Nov 25, 05
+  // Alice numbering scheme  - Jun 01,2006 
+  // IN:
+  // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
+  // nModule  - module number in SM,     0<= nModule  < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
+  // nIphi   - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv; 
+  // nIeta   - cell number in eta driection inside module; 0<= nIeta < fNETAdiv; 
+  // 
+ // OUT:
+  // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM
+  // ieta - have to change from 0 to (fNZ*fNETAdiv-1)
+  // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1)
+  //
+  static Int_t iphim, ietam;
+
+  GetModulePhiEtaIndexInSModule(nSupMod,nModule, iphim, ietam); 
+  //  ieta  = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM) 
+  ieta  = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM) 
+  iphi  = iphim*fNPHIdiv + nIphi;     // y(module) =  y(SM) 
+
+}
+
+
+//------------------------------------------------------------------------------------
+void AliJetDummyGeo::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule,  int &iphim, int &ietam) const
+{ 
+  // added nSupMod; - 19-oct-05 !
+  // Alice numbering scheme        - Jun 01,2006 
+  // ietam, iphi - indexes of module in two dimensional grid of SM
+  // ietam - have to change from 0 to fNZ-1
+  // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1)
+  static Int_t nphi;
+
+  if(nSupMod>=10) nphi = fNPhi/2;
+  else            nphi = fNPhi;
+
+  ietam = nModule/nphi;
+  iphim = nModule%nphi;
+}
+
+//------------------------------------------------------------------------------------
+Bool_t AliJetDummyGeo::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
+{
+  // Nov 17,2006
+  // stay here - phi problem as usual 
+  static Int_t nSupMod, i, ieta, iphi, etaShift, nphi;
+  static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc;
+  absId = nSupMod = - 1;
+  if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) {
+    // phi index first
+    phi    = TVector2::Phi_0_2pi(phi);
+    phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
+    nphi   = fPhiCentersOfCells.GetSize();
+    if(nSupMod>=10) {
+      phiLoc = phi - 190.*TMath::DegToRad();
+      nphi /= 2;
+    }
+
+    dmin   = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
+    iphi   = 0;
+    for(i=1; i<nphi; i++) {
+      d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
+      if(d < dmin) {
+        dmin = d;
+        iphi = i;
+      }
+    }
+    // odd SM are turned with respect of even SM - reverse indexes
+
+    // eta index
+    absEta   = TMath::Abs(eta);
+    etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
+    dmin     = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
+    ieta     = 0;
+    for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
+      d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
+      if(d < dmin) {
+        dmin = d;
+        ieta = i;
+      }
+    }
+
+    if(eta<0) iphi = (nphi-1) - iphi;
+    absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
+
+    return kTRUE;
+  }
+  return kFALSE;
+}
+
+//------------------------------------------------------------------------------------
+Bool_t AliJetDummyGeo::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
+{ 
+  // Return false if phi belongs a phi cracks between SM
+  static Int_t i;
+
+  if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
+
+  phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries
+  for(i=0; i<6; i++) {
+    if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) {
+      nSupMod = 2*i;
+      if(eta < 0.0) nSupMod++;
+      return kTRUE;
+    }
+  }
+  return kFALSE;
+}
+
+//------------------------------------------------------------------------------------
+Int_t  AliJetDummyGeo::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
+{
+  // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId
+  static Int_t ietam, iphim, nModule;
+  static Int_t nIeta, nIphi; // cell indexes in module
+
+  GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nModule);
+
+  nIeta = ieta%fNETAdiv;
+  nIeta = fNETAdiv - 1 - nIeta;
+  nIphi = iphi%fNPHIdiv;
+
+  return GetAbsCellId(nSupMod, nModule, nIphi, nIeta);
+}
+
+//------------------------------------------------------------------------------------
+void  AliJetDummyGeo::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta, 
+                       Int_t &iphim, Int_t &ietam, Int_t &nModule) const
+{
+  // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nModule)
+  static Int_t nphi;
+  nphi  = GetNumberOfModuleInPhiDirection(nSupMod);  
+
+  ietam  = ieta/fNETAdiv;
+  iphim  = iphi/fNPHIdiv;
+  nModule = ietam * nphi + iphim; 
+}
+
+//------------------------------------------------------------------------------------
+Int_t AliJetDummyGeo::GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const
+{ 
+  // 27-aug-04; 
+  // corr. 21-sep-04; 
+  //       13-oct-05; 110 degree case
+  // May 31, 2006; ALICE numbering scheme:
+  // 0 <= nSupMod < fNumberOfSuperModules
+  // 0 <= nModule  < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
+  // 0 <= nIphi   < fNPHIdiv
+  // 0 <= nIeta   < fNETAdiv
+  // 0 <= absid   < fNCells
+  static Int_t id=0; // have to change from 0 to fNCells-1
+  if(nSupMod >= 10) { // 110 degree case; last two supermodules
+    id  = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
+  } else {
+    id  = fNCellsInSupMod*nSupMod;
+  }
+  id += fNCellsInModule *nModule;
+  id += fNPHIdiv *nIphi;
+  id += nIeta;
+  if(id<0 || id >= fNCells) {
+    id = -TMath::Abs(id); // if negative something wrong
+  }
+  return id;
+}
+
+//------------------------------------------------------------------------------------
+Bool_t AliJetDummyGeo::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
+{
+  // 0<= nPhiSec <=4; phi in rad
+  // 0;  gap boundaries between  0th&2th  | 1th&3th SM
+  // 1;  gap boundaries between  2th&4th  | 3th&5th SM
+  // 2;  gap boundaries between  4th&6th  | 5th&7th SM
+  // 3;  gap boundaries between  6th&8th  | 7th&9th SM
+  // 4;  gap boundaries between  8th&10th | 9th&11th SM
+  if(nPhiSec<0 || nPhiSec >4) return kFALSE; 
+  phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
+  phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
+  return kTRUE; 
+}
+
+//------------------------------------------------------------------------------------
+void  AliJetDummyGeo::GetTransformationForSM()
+{
+  // Uses the geometry manager to load the transformation matrix
+  // for the supermodules
+  // Unused after 19 Jan, 2007 - keep for compatibility; 
+
+  return;
+  static Bool_t transInit=kFALSE;
+  if(transInit) return;
+
+  int i=0;
+  if(gGeoManager == 0) {
+    Info("CreateTransformationForSM() "," Load geometry : TGeoManager::Import()");
+    assert(0);
+  }
+
+  TGeoNode *tn = gGeoManager->GetTopNode();
+  TGeoNode *node=0, *xen1 = 0;
+  for(i=0; i<tn->GetNdaughters(); i++) {
+    node = tn->GetDaughter(i);
+    TString ns(node->GetName());
+    if(ns.Contains(GetNameOfEMCALEnvelope())) {
+      xen1 = node;
+      break;
+    }
+  }
+
+  if(!xen1) {
+    Info("CreateTransformationForSM() "," geometry has not EMCAL envelope with name %s", 
+    GetNameOfEMCALEnvelope());
+    assert(0);
+  }
+  printf(" i %i : EMCAL Envelope is %s : #SM %i \n", i, xen1->GetName(), xen1->GetNdaughters());
+  for(i=0; i<xen1->GetNdaughters(); i++) {
+    TGeoNodeMatrix *sm = (TGeoNodeMatrix*)xen1->GetDaughter(i);
+    fMatrixOfSM[i] = sm->GetMatrix();
+  }
+  printf("transInit %d: ", transInit);
+  transInit = kTRUE;
+}
+
index 475d5ad..2d29bdd 100644 (file)
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
  
+//
+// Temporarily added to define part of the EMCal geometry
+// necessary for the jet finder
+//
+
 #include <TObject.h>
+#include <TArrayD.h>
+#include <TMath.h>
+#include <TVector3.h>
+
+class TGeoMatrix;
 
 class AliJetDummyGeo : public TObject 
 {
  public: 
-    AliJetDummyGeo(){;}
-    virtual ~AliJetDummyGeo(){;}
-    static AliJetDummyGeo* GetInstance() {return new AliJetDummyGeo();}
-    static AliJetDummyGeo* GetInstance(char* /*name*/, char* /*title*/)
-       {return new AliJetDummyGeo();}
-    Int_t GetNCells(){ return 0;}
-    Float_t GetArm1EtaMin() {return 0.;}
-    Float_t GetArm1EtaMax() {return 0.;}
-    Float_t GetArm1PhiMin() {return 0.;}
-    Float_t GetArm1PhiMax() {return 0.;}
-    void    EtaPhiFromIndex(Int_t /*id*/, Float_t& /*eta*/, Float_t& /*phi*/)
-       {;}
-    Int_t   TowerIndexFromEtaPhi2(Float_t /*eta*/, Float_t /*phi*/) {return 0;}
-    void    GetTransformationForSM(){;}
-    Float_t GetSampling() {return 0.;}
-    ClassDef(AliJetDummyGeo,1)
+  AliJetDummyGeo();
+  AliJetDummyGeo(const AliJetDummyGeo& geom);
+  virtual ~AliJetDummyGeo();
+  static AliJetDummyGeo* GetInstance() {return new AliJetDummyGeo();}
+  static AliJetDummyGeo* GetInstance(char* /*name*/, char* /*title*/)
+    {return new AliJetDummyGeo();}
+  Char_t* GetNameOfEMCALEnvelope() const {return "XEN1";}
+  Float_t GetEnvelop(Int_t index) const { return fEnvelop[index];}
+  Float_t AngleFromEta(Float_t eta) const { 
+    // returns theta in radians for a given pseudorapidity
+    return 2.0*TMath::ATan(TMath::Exp(-eta));
+  }
+  Float_t ZFromEtaR(Float_t r,Float_t eta) const { 
+    // returns z in for a given
+    // pseudorapidity and r=sqrt(x*x+y*y).
+    return r/TMath::Tan(AngleFromEta(eta));
+  }
+  Int_t   GetNCells()     {return fNCells;}
+  Float_t GetPhiModuleSize() const  {return fPhiModuleSize;}
+  Float_t GetEtaModuleSize() const  {return fEtaModuleSize;}
+  Float_t GetShellThickness() const {return fShellThickness;}
+  Int_t   GetNPhi() const {return fNPhi;}
+  Int_t   GetNumberOfSuperModules() const {return fNumberOfSuperModules;}
+  Float_t GetArm1EtaMin() {return fArm1EtaMin;}
+  Float_t GetArm1EtaMax() {return fArm1EtaMax;}
+  Float_t GetArm1PhiMin() {return fArm1PhiMin;}
+  Float_t GetArm1PhiMax() {return fArm1PhiMax;}
+  void    EtaPhiFromIndex(Int_t /*id*/, Float_t& /*eta*/, Float_t& /*phi*/);
+  void    GetGlobal(const Double_t *loc, Double_t *glob, int ind) const;
+  void    GetGlobal(Int_t absId, Double_t glob[3]) const;
+  void    GetGlobal(Int_t absId, TVector3 &vglob) const;
+  Bool_t  RelPosCellInSModule(Int_t absId, Double_t loc[3]) const;
+  Bool_t  RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const;
+  Bool_t  CheckAbsCellId(Int_t absId) const;
+  Bool_t  GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const;
+  void    GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta, int &iphi, int &ieta) const;
+  void    GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule,  int &iphim, int &ietam) const;
+  Bool_t  GetAbsCellIdFromEtaPhi(Double_t eta,Double_t phi, Int_t &absId) const;
+  Bool_t  SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const;
+  Int_t   GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const;
+  void    GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta, 
+                                                  Int_t &iphim, Int_t &ietam, Int_t &nModule) const;
+  Int_t   GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const;
+  Int_t   GetNumberOfModuleInPhiDirection(Int_t nSupMod)  const
+  {
+    // inline function
+    if(nSupMod>=10) return fNPhi/2;
+    else            return fNPhi;
+  }
+  Bool_t  GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const;
+  void    GetTransformationForSM();
+  Float_t GetSampling() {return fSampling;}
+
+ protected:
+  Float_t fArm1PhiMin;                         // Minimum angular position of EMCAL in Phi (degrees)
+  Float_t fArm1PhiMax;                 // Maximum angular position of EMCAL in Phi (degrees)
+  Float_t fArm1EtaMin;                 // Minimum pseudorapidity position of EMCAL in Eta
+  Float_t fArm1EtaMax;                         // Maximum pseudorapidity position of EMCAL in Eta
+  Int_t   fNumberOfSuperModules;
+  Float_t fSteelFrontThick;             // Thickness of the front stell face of the support box - 9-sep-04
+  Float_t fEnvelop[3];                 // the GEANT TUB for the detector 
+  Float_t fIPDistance;                  // Radial Distance of the inner surface of the EMCAL
+  Float_t fZLength;                    // Total length in z direction
+  Float_t fPhiGapForSM;                 // Gap betweeen supermodules in phi direction
+  Int_t   fNPhi;                       // Number of Towers in the PHI direction  
+  Int_t   fNZ;                         // Number of Towers in the Z direction
+  Float_t fPhiModuleSize;               // Phi -> X
+  Float_t fEtaModuleSize;               // Eta -> Y
+  Int_t   fNPHIdiv;                     // number phi divizion of module
+  Int_t   fNETAdiv;                     // number eta divizion of module
+  Int_t   fNECLayers;                   // number of scintillator layers
+  Float_t fECScintThick;                // cm, Thickness of the scintillators
+  Float_t fECPbRadThickness;            // cm, Thickness of the Pb radiators
+  Float_t fSampling;                    // Sampling factor
+  Float_t fTrd1Angle;                   // angle in x-z plane (in degree) 
+  Int_t   fNCellsInModule;              // number cell in module
+  Int_t   fNCellsInSupMod;              // number cell in super module
+  Int_t   fNCells;                      // number of cells in calo
+  Float_t fLongModuleSize;              // Size of long module
+  Float_t f2Trd1Dx2;                    // 2*dx2 for TRD1
+  Float_t fShellThickness;             // Total thickness in (x,y) direction
+  Float_t fEtaMaxOfTRD1;                // max eta in case of TRD1 geometry (see AliEMCALShishKebabTrd1Module)
+  Float_t fParSM[3];                    // SM sizes as in GEANT (TRD1)
+  TArrayD fPhiBoundariesOfSM;           // phi boundaries of SM in rad; size is fNumberOfSuperModules;
+  TArrayD fPhiCentersOfSM;              // phi of centers of SMl size is fNumberOfSuperModules/2
+  TGeoMatrix* fMatrixOfSM[12];          //![fNumberOfSuperModules]; get from gGeoManager;
+  TArrayD fCentersOfCellsEtaDir;        // size fNEta*fNETAdiv (for TRD1 only) (eta or z in SM, in cm)
+  TArrayD fCentersOfCellsXDir;          // size fNEta*fNETAdiv (for TRD1 only) (       x in SM, in cm)
+  TArrayD fCentersOfCellsPhiDir;        // size fNPhi*fNPHIdiv (for TRD1 only) (phi or y in SM, in cm)
+  TArrayD fEtaCentersOfCells;           // [fNEta*fNETAdiv*fNPhi*fNPHIdiv], positive direction (eta>0); 
+                                        // eta depend from phi position; 
+  TArrayD fPhiCentersOfCells;           // [fNPhi*fNPHIdiv] from center of SM (-10. < phi < +10.)
+
+  ClassDef(AliJetDummyGeo,1)
 };
  
 #endif
index 3975f25..733b26a 100755 (executable)
 //          Magali Estienne <magali.estienne@IReS.in2p3.fr>
 //------------------------------------------------------------------------- 
 
-
+// --- Standard library ---
 #include <Riostream.h>
+
+// --- ROOT system ---
 #include <TSystem.h>
+#include <TStopwatch.h>
 #include <TLorentzVector.h>
 #include <TVector3.h>
+#include <TTask.h>
 #include <TGeoManager.h>
+#include <assert.h>
+#include <TRefArray.h>
+#include <TMath.h>
+#include <TChain.h>
 
+
+// --- AliRoot header files ---
 #include "AliJetESDReader.h"
 #include "AliJetESDReaderHeader.h"
 #include "AliESDEvent.h"
 #include "AliESD.h"
 #include "AliESDtrack.h"
-//#include "AliEMCALGeometry.h"
 #include "AliJetDummyGeo.h"
 #include "AliJetFillUnitArrayTracks.h"
+#include "AliJetFillUnitArrayEMCalDigits.h"
 #include "AliJetUnitArray.h"
 
 ClassImp(AliJetESDReader)
 
 AliJetESDReader::AliJetESDReader():
-    AliJetReader(),  
-    fGeom(0),
-    fChain(0x0),
-    fESD(0x0),
-    fHadCorr(0x0),
-    fTpcGrid(0x0),
-    fEmcalGrid(0x0),
-    fPtCut(0),
-    fHCorrection(0),
-    fNumUnits(0),
-    fDebug(0),
-    fNIn(0),
-    fOpt(0),
-    fNeta(0),
-    fNphi(0),
-    fArrayInitialised(0) 
+  AliJetReader(),  
+  fGeom(0),
+  fChain(0x0),
+  fESD(0x0),
+  fHadCorr(0x0),
+  fTpcGrid(0x0),
+  fEmcalGrid(0x0),
+  fGrid0(0),
+  fGrid1(0),
+  fGrid2(0),
+  fGrid3(0),
+  fGrid4(0),
+  fPtCut(0),
+  fHCorrection(0),
+  fNumUnits(0),
+  fDebug(0),
+  fMass(0),
+  fSign(0),
+  fNIn(0),
+  fOpt(0),
+  fDZ(0),
+  fNeta(0),
+  fNphi(0),
+  fArrayInitialised(0)
 {
-  // Constructor    
+  // Constructor 
 }
 
 //____________________________________________________________________________
-
 AliJetESDReader::~AliJetESDReader()
 {
   // Destructor
@@ -69,10 +86,17 @@ AliJetESDReader::~AliJetESDReader()
     delete fESD;
     delete fTpcGrid;
     delete fEmcalGrid;
+    if(fDZ)
+      {
+       delete fGrid0;
+       delete fGrid1;
+       delete fGrid2;
+       delete fGrid3;
+       delete fGrid4;
+      }
 }
 
 //____________________________________________________________________________
-
 void AliJetESDReader::OpenInputFiles()
 {
   // chain for the ESDs
@@ -82,8 +106,7 @@ void AliJetESDReader::OpenInputFiles()
    const char* dirName=fReaderHeader->GetDirectory();
    const char* pattern=fReaderHeader->GetPattern();
     
-//   // Add files matching patters to the chain
-  
+   // Add files matching patters to the chain
    void *dir  = gSystem->OpenDirectory(dirName);
    const char *name = 0x0;
    int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
@@ -92,19 +115,20 @@ void AliJetESDReader::OpenInputFiles()
        if (a>=nesd) continue;
        
        if (strstr(name,pattern)){
-          char path[256];
-          sprintf(path,"%s/%s/AliESDs.root",dirName,name);
-          fChain->AddFile(path);
-          a++;
+        printf("Adding %s\n",name);
+        char path[256];
+        //        sprintf(path,"%s/%s/AliESDs.root",dirName,name);
+        sprintf(path,"%s/%s/AliESDs.root",dirName,name);
+        fChain->AddFile(path);
+        a++;
        }
    }
   
   gSystem->FreeDirectory(dir);
   
+  fESD = new AliESDEvent();
+  fESD->ReadFromTree(fChain);
 
-  fESD = 0;
-  fChain->SetBranchAddress("ESD",    &fESD);
-  
   int nMax = fChain->GetEntries(); 
 
   printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
@@ -118,13 +142,13 @@ void AliJetESDReader::OpenInputFiles()
   }
 }
 
+//____________________________________________________________________________
 void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) {
     // Connect the tree
      fESD   = (AliESDEvent*) esd;
 }
 
 //____________________________________________________________________________
-
 Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/)
 {
   // Fill momentum array
@@ -137,11 +161,11 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/)
   // clear momentum array
   ClearArray();
   fDebug = fReaderHeader->GetDebug();
-  
+  fOpt = fReaderHeader->GetDetector();
+
   if (!fESD) {
       return kFALSE;
   }
-  
 
   // get number of tracks in event (for the loop)
   nt = fESD->GetNumberOfTracks();
@@ -184,50 +208,138 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/)
   // set the signal flags
   fSignalFlag.Set(goodTrack,sflag);
   fCutFlag.Set(goodTrack,cflag);
+  return kTRUE;
 
-//
-//
-  if (fTpcGrid || fEmcalGrid) {
-      SetEMCALGeometry();
-      InitParameters();
-      AliJetFillUnitArrayTracks *fillUAFromTracks = new AliJetFillUnitArrayTracks(); 
-      fillUAFromTracks->SetReaderHeader(fReaderHeader);
-      fillUAFromTracks->SetMomentumArray(fMomentumArray);
-      fillUAFromTracks->SetTPCGrid(fTpcGrid);
-      fillUAFromTracks->SetEMCalGrid(fEmcalGrid);
-      fillUAFromTracks->SetHadCorrection(fHCorrection);
-      fillUAFromTracks->SetHadCorrector(fHadCorr);
-      fNeta = fillUAFromTracks->GetNeta();
-      fNphi = fillUAFromTracks->GetNphi();
-      fillUAFromTracks->SetActive(kFALSE);
-      // TPC only or Digits+TPC or Clusters+TPC
-      if(fOpt%2==!0 && fOpt!=0) { 
-         fillUAFromTracks->SetActive(kTRUE);
-         fillUAFromTracks->SetUnitArray(fUnitArray);
-         fillUAFromTracks->ExecuteTask("tpc");
-      }
+}
+
+//____________________________________________________________________________
+void AliJetESDReader::CreateTasks()
+{
+  fDebug = fReaderHeader->GetDebug();
+  fDZ = fReaderHeader->GetDZ();
+
+  // Init EMCAL geometry and create UnitArray object
+  SetEMCALGeometry();
+  InitParameters();
+  InitUnitArray();
+
+  fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
+  fFillUAFromTracks = new AliJetFillUnitArrayTracks(); 
+  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->SetHadCorrection(fHCorrection);
+  fFillUAFromTracks->SetHadCorrector(fHadCorr);
+  fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits();
+  fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
+  fFillUAFromEMCalDigits->SetGeom(fGeom);
+  fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
+  fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
+  //      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;
+   
+}
+
+//____________________________________________________________________________
+//void AliJetESDReader::ExecTasks(Int_t event)
+Bool_t AliJetESDReader::ExecTasks(Int_t /*event*/)
+{
+  // clear momentum array
+  Int_t nEntRef = fRefArray->GetEntries();
+
+  for(Int_t i=0; i<nEntRef; i++)
+    {      
+      ((AliJetUnitArray*)fRefArray->At(i))->SetUnitTrackID(0);
+      ((AliJetUnitArray*)fRefArray->At(i))->SetUnitEnergy(0.);
+      ((AliJetUnitArray*)fRefArray->At(i))->SetUnitCutFlag(kPtSmaller);
+      ((AliJetUnitArray*)fRefArray->At(i))->SetUnitDetectorFlag(kTpc);
+      ((AliJetUnitArray*)fRefArray->At(i))->SetUnitFlag(kOutJet);
+    }
+
+  ClearArray();
+
+  fDebug = fReaderHeader->GetDebug();
+  fOpt = fReaderHeader->GetDetector();
+  //  InitParameters();
+
+  if(!fESD) {
+    return kFALSE;
+  }
   
-      delete fillUAFromTracks;
+  /*
+  // get event from chain
+  // For TSelectors
+  //  fChain->GetTree()->GetEntry(event);
+  // For interactive process
+  //  fChain->GetEntry(event);
+  fChain->GetEvent(event);
+  */
+
+  // TPC only or Digits+TPC or Clusters+TPC
+  if(fOpt%2==!0 && fOpt!=0){ 
+    fFillUAFromTracks->SetESD(fESD);
+    fFillUAFromTracks->SetActive(kTRUE);
+    fFillUAFromTracks->SetUnitArray(fUnitArray);
+    fFillUAFromTracks->SetRefArray(fRefArray);
+    //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed  !!!
+                                                 // Incompatibility with Andrei's analysis framework
+    fFillUAFromTracks->Exec("tpc");
+    if(fOpt==1){
+      fNumCandidate = fFillUAFromTracks->GetMult();
+      fNumCandidateCut = fFillUAFromTracks->GetMultCut();
+    }
+  }
+
+  // Digits only or Digits+TPC  
+  if(fOpt>=2 && fOpt<=3){
+    fFillUAFromEMCalDigits->SetESD(fESD);
+    fFillUAFromEMCalDigits->SetActive(kTRUE);
+    fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
+    fFillUAFromEMCalDigits->SetRefArray(fRefArray);
+    fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
+    fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
+    fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily changed !!!
+    fNumCandidate = fFillUAFromEMCalDigits->GetMult();
+    fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
   }
 
+  //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
+
   return kTRUE;
 }
 
-
+//____________________________________________________________________________
 void AliJetESDReader::SetEMCALGeometry()
 {
   // Define EMCAL geometry to be able to read ESDs
     fGeom = AliJetDummyGeo::GetInstance();
     if (fGeom == 0)
        fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
-    
+
     // To be setted to run some AliEMCALGeometry functions
     TGeoManager::Import("geometry.root");
-    fGeom->GetTransformationForSM();  
+    //    fGeom->GetTransformationForSM();  
     printf("\n EMCal Geometry set ! \n");
 
 }
-  
+
+
+//____________________________________________________________________________  
 void AliJetESDReader::InitParameters()
 {
     // Initialise parameters
@@ -237,16 +349,153 @@ void AliJetESDReader::InitParameters()
     if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 }
 
+//____________________________________________________________________________
 void AliJetESDReader::InitUnitArray()
 {
   //Initialises unit arrays
-    Int_t nElements = fTpcGrid->GetNEntries();
-    if(fArrayInitialised) delete [] fUnitArray;
-    if(fTpcGrid->GetGridType()==0)
-       fUnitArray = new AliJetUnitArray[nElements];
-    if(fTpcGrid->GetGridType()==1)
-       fUnitArray = new AliJetUnitArray[fNumUnits+nElements];
-    fArrayInitialised = 1;
+  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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
+       }
+    }
+
+  if(fTpcGrid->GetGridType()==1)
+    {
+      Int_t nGaps = 0;
+      Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
+
+      if(fDZ)
+       {
+         // Define a grid of cell for the gaps between SM
+         Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
+         Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
+         fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
+         fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
+         fGrid0->SetGridType(0);
+         fGrid0->SetMatrixIndexes();
+         fGrid0->SetIndexIJ();
+         n0 = fGrid0->GetNEntries();
+         fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
+         fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
+         fGrid1->SetGridType(0);
+         fGrid1->SetMatrixIndexes();
+         fGrid1->SetIndexIJ();
+         n1 = fGrid1->GetNEntries();
+         fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
+         fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
+         fGrid2->SetGridType(0);
+         fGrid2->SetMatrixIndexes();
+         fGrid2->SetIndexIJ();
+         n2 = fGrid2->GetNEntries();
+         fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
+         fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
+         fGrid3->SetGridType(0);  
+         fGrid3->SetMatrixIndexes();
+         fGrid3->SetIndexIJ();
+         n3 = fGrid3->GetNEntries();
+         fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
+         fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
+         fGrid4->SetGridType(0);
+         fGrid4->SetMatrixIndexes();
+         fGrid4->SetIndexIJ();
+         n4 = fGrid4->GetNEntries();
+
+         if(fDebug>1) 
+           {
+             cout << "n0 cells: " << n0 << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl;
+             cout << "n1 cells: " << n1 << "phimin1: " << phimin1 << ", phimax1: " << phimax1 << endl;
+             cout << "n2 cells: " << n2 << "phimin2: " << phimin2 << ", phimax2: " << phimax2 << endl;
+             cout << "n3 cells: " << n3 << "phimin3: " << phimin3 << ", phimax3: " << phimax3 << endl;
+             cout << "n4 cells: " << n4 << "phimin4: " << phimin4 << ", phimax4: " << phimax4 << endl;
+           }
+         
+         nGaps = n0+n1+n2+n3+n4;
+
+       }
+
+      for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++) 
+       {
+         if(nBin<fNumUnits)
+           {
+             fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry 
+             // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
+             phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
+             Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
+             Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
+             new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,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.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
+                   }
+               }
+             } // end if(fDZ)
+           } // end else 2
+         } // end else 1
+       } // end loop on nBin
+    } // end grid type == 1
+  fArrayInitialised = 1;
 }
 
 
+
index d5cdf59..55d41cb 100755 (executable)
@@ -18,7 +18,6 @@
 #include "AliJetUnitArray.h"
 #include "AliJetGrid.h"
 class AliJetESDReaderHeader;
-class AliEMCALGeometry;
 class AliJetDummyGeo;
 class AliJetHadronCorrection;
 class AliJetUnitArray;
@@ -30,10 +29,18 @@ class AliJetESDReader : public AliJetReader
  public: 
   AliJetESDReader();
   virtual ~AliJetESDReader();
+
+  // Getters
+  Float_t GetTrackMass() const {return fMass;}  // returns mass of the track
+  Int_t   GetTrackSign() const {return fSign;}  // returns sign of the track
+
   // Setters
   Bool_t FillMomentumArray(Int_t event); 
   void   OpenInputFiles();
   void   InitUnitArray();
+  void   CreateTasks();
+  //  void   ExecTasks(Int_t event);
+  Bool_t   ExecTasks(Int_t event);
   void   SetInputEvent(TObject* esd, TObject* aod, TObject* mc);
   virtual void SetTPCGrid(AliJetGrid *grid)   {fTpcGrid = grid;}
   virtual void SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;}
@@ -50,19 +57,25 @@ class AliJetESDReader : public AliJetReader
   AliJetHadronCorrectionv1   *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                       fNumUnits;         // Number of units in the unit object array
                                                  // (same as num towers in EMCAL)
   Int_t                       fDebug;            // Debug option
+  Float_t                     fMass;    // Particle mass
+  Int_t                       fSign;    // Particle sign
   Int_t                       fNIn;              // Number of Array filled in UnitArray
   Int_t                       fOpt;              // Detector to be used for jet reconstruction
+  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
   
-
-
   ClassDef(AliJetESDReader,1)
 };
  
diff --git a/JETAN/AliJetFillUnitArrayEMCalDigits.cxx b/JETAN/AliJetFillUnitArrayEMCalDigits.cxx
new file mode 100644 (file)
index 0000000..a078364
--- /dev/null
@@ -0,0 +1,248 @@
+
+/**************************************************************************
+ * 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 with EMCal information 
+// Called by ESD reader for jet analysis
+// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
+//------------------------------------------------------------- 
+
+// --- Standard library ---
+#include <Riostream.h>
+
+// --- ROOT system ---
+#include <TSystem.h>
+#include <TLorentzVector.h>
+#include <TRefArray.h> 
+#include "TTask.h"
+#include <TGeoManager.h>
+#include <TMath.h>
+#include <TClonesArray.h>
+#include <TArrayS.h>
+
+// --- AliRoot header files ---
+#include "AliJetFinder.h"
+#include "AliJetReaderHeader.h"
+#include "AliJetReader.h"
+#include "AliJetESDReader.h"
+#include "AliJetESDReaderHeader.h"
+//#include "AliESD.h"
+#include "AliESDEvent.h"
+#include "AliJetDummyGeo.h"
+#include "AliESDCaloCluster.h"
+#include "AliJetUnitArray.h"
+#include "AliJetFillUnitArrayEMCalDigits.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliCDBEntry.h"
+
+ClassImp(AliJetFillUnitArrayEMCalDigits)
+
+//_____________________________________________________________________________
+AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits()
+  : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
+    fESD(0),    
+    fNumUnits(0), 
+    fEtaMinCal(0.),
+    fEtaMaxCal(0.),
+    fPhiMinCal(0.),
+    fPhiMaxCal(0.),
+    fNIn(0),
+    fOpt(0),
+    fDebug(0),
+    fNCEMCAL(0),
+    fNCPHOS(0),
+    fNCCalo(0),
+    fTPCGrid(0x0),
+    fEMCalGrid(0x0),
+    fReaderHeader(0x0),
+    fMomentumArray(0x0),
+    fUnitArray(0x0),
+    fRefArray(0x0),
+    fGeom(0x0),
+    fClus(0x0),
+    fNDigitEmcal(0),
+    fNDigitEmcalCut(0)
+{
+  // constructor
+}
+
+//_____________________________________________________________________________
+AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent */*esd*/)
+  : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
+    fESD(0),    
+    fNumUnits(0), 
+    fEtaMinCal(0.),
+    fEtaMaxCal(0.),
+    fPhiMinCal(0.),
+    fPhiMaxCal(0.),
+    fNIn(0),
+    fOpt(0),
+    fDebug(0),
+    fNCEMCAL(0),
+    fNCPHOS(0),
+    fNCCalo(0),
+    fTPCGrid(0x0),
+    fEMCalGrid(0x0),
+    fReaderHeader(0x0),
+    fMomentumArray(0x0),
+    fUnitArray(0x0),
+    fRefArray(0x0),
+    fGeom(0x0),
+    fClus(0x0),
+    fNDigitEmcal(0),
+    fNDigitEmcalCut(0)
+{
+  // constructor
+}
+
+
+//____________________________________________________________________________
+void AliJetFillUnitArrayEMCalDigits::InitParameters()
+{
+  fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
+
+  fEtaMinCal = fGeom->GetArm1EtaMin();
+  fEtaMaxCal = fGeom->GetArm1EtaMax();
+  fPhiMinCal = fGeom->GetArm1PhiMin();
+  fPhiMaxCal = fGeom->GetArm1PhiMax(); 
+  fClus      = 0;
+
+  if(fDebug>1) printf("\n EMCAL parameters initiated ! \n");
+
+}
+
+//_____________________________________________________________________________
+AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits()
+{
+  // destructor
+}
+
+//_____________________________________________________________________________
+void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/)
+{
+  //
+  // Main method.
+  // Explain
+
+  fDebug = fReaderHeader->GetDebug();
+  fOpt = fReaderHeader->GetDetector();
+
+  // Init parameters
+  InitParameters();
+
+  // Get number of clusters from EMCAL
+  
+  Int_t   nDigitTot      = 0;
+  Int_t   goodDigit      = 0;
+  Int_t   beg            = 0;
+  Int_t   end            = 0;
+  Float_t ptMin          = fReaderHeader->GetPtCut();
+
+  // Loop over calo clusters 
+  //------------------------------------------------------------------
+  Int_t type = 0; 
+  Int_t index = 0;
+
+  // Total number of EMCAL cluster
+  end =  fESD->GetNumberOfCaloClusters();
+
+  for(Int_t j = beg; j < end; j++) {
+      fClus = fESD->GetCaloCluster(j);
+      if(!fClus->IsEMCAL()) continue;
+
+      type = fClus->GetClusterType(); 
+      index = fClus->GetID();
+      nDigitTot = fClus->GetNumberOfDigits();
+      
+      // Keep clusters or pseudo clusters
+      if (type != AliESDCaloCluster::kClusterv1) continue;
+      //      if (type != AliESDCaloCluster::kPseudoCluster) continue;
+
+      // Get the digit index and the digit information
+      //============================================================
+         
+      // Get number of digits in a cluster
+      Int_t nD = fClus->GetNumberOfDigits();
+
+      TArrayS *digID = fClus->GetDigitIndex();      
+      TArrayS *digEnergy = fClus->GetDigitAmplitude();
+      Float_t *digitEnergy = new Float_t[nD];      
+      //      Float_t digitEn = 0.;
+         
+      // Loop over digits
+      for(Int_t k=0; k<nD; k++) {
+       
+       // Convert energy in GeV
+       Int_t idF = (Int_t)digID->At(k);
+       // Calibration for an energy in GeV
+       digitEnergy[k] = (Float_t)digEnergy->At(k)/500.; 
+
+       // Second method to extract eta, phi positions of a digit
+       //=================================================================
+
+       Float_t etaD=-10., phiD=-10.;
+       fGeom->EtaPhiFromIndex(idF,etaD,phiD); 
+       //          fEMCalGrid->GetEtaPhiFromIndex2(idF,phiD,etaD);
+       phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
+
+       Float_t etDigit = digitEnergy[k]*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
+
+       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idF);
+       if(uArray->GetUnitEnergy() == 0.) goodDigit++;
+
+       Float_t unitEnergy = 0.;
+       Bool_t ok = kFALSE;
+       unitEnergy = uArray->GetUnitEnergy();
+       if(unitEnergy==0){
+         fRefArray->Add(uArray);
+         fNDigitEmcal++;
+         ok = kTRUE;
+       }
+       uArray->SetUnitEnergy(unitEnergy+etDigit);
+       // Put a cut flag
+       if(uArray->GetUnitEnergy()<ptMin)
+         uArray->SetUnitCutFlag(kPtSmaller);
+       else {
+         uArray->SetUnitCutFlag(kPtHigher);
+         if(ok) fNDigitEmcalCut++;
+       }
+       // Detector flag
+       if(unitEnergy>0.) 
+         uArray->SetUnitDetectorFlag(kAll);
+       else uArray->SetUnitDetectorFlag(kEmcal);
+       
+       // This is for jet multiplicity
+       uArray->SetUnitClusterID(index);
+       
+       if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
+       
+      } // End loop over digits
+      
+  } // End loop over clusters
+  
+  fNIn += goodDigit;
+
+}
+
+//_____________________________________________________________________________
+Float_t  AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg)
+{
+  //  return (180./TMath::Pi())*2.*atan(exp(-arg));
+  return 2.*atan(exp(-arg));
+
+
+}
diff --git a/JETAN/AliJetFillUnitArrayEMCalDigits.h b/JETAN/AliJetFillUnitArrayEMCalDigits.h
new file mode 100644 (file)
index 0000000..e59332c
--- /dev/null
@@ -0,0 +1,95 @@
+#ifndef ALIJETFILLUNITARRAYEMCALDIGITS_H
+#define ALIJETFILLUNITARRAYEMCALDIGITS_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)
+//---------------------------------------------------------------------
+
+#ifndef ROOT_TTask
+#include "TTask.h"
+#endif
+
+class AliJetDummyGeo;
+class AliESDCaloCluster;
+class AliEMCALCalibData;
+class AliJetReader;
+class AliJetESDReader;
+class TClonesArray;
+class TRefArray;
+class AliJetUnitArray;
+//class AliESD;
+class AliESDEvent;
+class AliJetGrid;
+
+class AliJetFillUnitArrayEMCalDigits : public TTask
+{
+ public: 
+  AliJetFillUnitArrayEMCalDigits();
+  AliJetFillUnitArrayEMCalDigits(Int_t event);
+  AliJetFillUnitArrayEMCalDigits(AliESD *fESD);
+  AliJetFillUnitArrayEMCalDigits(AliESDEvent *fESD);
+  virtual ~AliJetFillUnitArrayEMCalDigits();
+  
+  // Setter
+  void SetReaderHeader(AliJetReaderHeader *readerHeader) {fReaderHeader = readerHeader;}
+  void SetGeom(AliJetDummyGeo *geom) {fGeom = geom;}
+  void SetMomentumArray(TClonesArray *momentumArray) {fMomentumArray = momentumArray;}
+  void SetUnitArray(TClonesArray *unitArray) {fUnitArray = unitArray;}
+  void SetRefArray(TRefArray *refArray) {fRefArray = refArray;}
+  void SetTPCGrid(AliJetGrid *grid) {fTPCGrid = grid;}
+  void SetEMCalGrid(AliJetGrid *grid) {fEMCalGrid = grid;}
+  //  void SetESD(AliESD *esd) {fESD = esd;}
+  void SetESD(AliESDEvent *esd) {fESD = esd;}
+  void SetInitMult(Int_t mult) {fNDigitEmcal = mult;}
+  void SetInitMultCut(Int_t multcut) {fNDigitEmcalCut = multcut;}
+
+  // Getter
+  TClonesArray* GetUnitArray() {return fUnitArray;}
+  TRefArray*    GetRefArray()  {return fRefArray;}
+  Int_t         GetMult()      {return fNDigitEmcal;}
+  Int_t         GetMultCut()   {return fNDigitEmcalCut;}
+
+  // Other
+  void          Exec(Option_t*);
+  Float_t       EtaToTheta(Float_t arg);
+ private:
+  void InitParameters();
+
+ protected:
+  AliESDEvent  *fESD; // ESD
+  Int_t   fNumUnits;  // Number of units in the unit object array (same as num towers in EMCAL)
+  Float_t fEtaMinCal; // Define EMCAL acceptance in Eta
+  Float_t fEtaMaxCal; // Define EMCAL acceptance in Eta
+  Float_t fPhiMinCal; // Define EMCAL acceptance in Phi
+  Float_t fPhiMaxCal; // Define EMCAL acceptance in Phi
+  Int_t   fNIn;       // Number of Array filled in UnitArray
+  Int_t   fOpt;       // Detector to be used for jet reconstruction
+  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
+
+  AliJetGrid         *fTPCGrid;           // Define filled grid
+  AliJetGrid         *fEMCalGrid;         // Define filled grid
+
+  AliJetReaderHeader *fReaderHeader;      // ReaderHeader
+  TClonesArray       *fMomentumArray;     // MomentumArray
+  TClonesArray       *fUnitArray;         // UnitArray
+  TRefArray          *fRefArray;          // UnitArray
+  AliJetDummyGeo     *fGeom;              // Set EMCal geometry
+  
+  AliESDCaloCluster *fClus;               //! 
+  Int_t fNDigitEmcal;                     //!
+  Int_t fNDigitEmcalCut;                  //!
+
+  
+
+  ClassDef(AliJetFillUnitArrayEMCalDigits,1) // Fill Unit Array with tpc and/or emcal information
+};
+
+#endif
index 1cf9d0e..e5cdc2e 100644 (file)
@@ -1,4 +1,3 @@
-
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-//======================================================================
-// *** November 2006
-// Author: magali.estienne@ires.in2p3.fr
-// 1) Define 2 grids and take (eta,phi) from the grid or use the grid for the TPC and  
-// EtaPhiFromIndex and TowerIndexFromEtaPhi for the particles in EMCAL acceptance
-//     2 options are possible : fGrid==0, work only with full TPC acceptance (for the moment)
-//                              fGrid==1, work with a part of the TPC acceptance  
-// 2) Try to implement 2 full grids for TPC and EMCal separately and to merge them
-// 3) Need to include Dead-zone -> Wait for exact positions in the new detector geometry
-// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
-//======================================================================
-// ***September 2006
-// TTask : Fill Unit Array for the Tracks information 
-// Called by ESD reader for jet analysis
-// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
-//======================================================================
-// *** July 2006
-// 1) When the tracks are in the EMCal acceptance, the functions EtaPhiFromIndex
-// and TowerIndexFromEtaPhi in the AliEMCALGeometry class are used to extract the
-// index or the eta, phi position of a grid.
-// 2) Define a grid for TPC
-// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
-//======================================================================
-// ***July 2006
+
+//----------------------------------------------------------------------
 // Fill Unit Array class 
-// Class used by AliJetESDReader to fill a UnitArray from the information extracted 
-// from the particle tracks
+// Class used by AliJetESDReader to fill a UnitArray from the information
+// extracted from the particle tracks
 // Author: magali.estienne@ires.in2p3.fr
-//======================================================================
+//----------------------------------------------------------------------
 
 
 // --- Standard library ---
 // --- ROOT system ---
 #include <TSystem.h>
 #include <TLorentzVector.h>
+#include <TRefArray.h> 
 #include <TVector3.h>
-//#include "Math/Vector3D.h"
-//#include "Math/Vector3Dfwd.h"
 #include "TTask.h"
 #include <TGeoManager.h>
 #include <TMatrixD.h>
 #include <TArrayD.h>
+#include <TMath.h>
+#include <TClonesArray.h>
 
 // --- AliRoot header files ---
 #include "AliJetFinder.h"
@@ -65,8 +43,8 @@
 #include "AliJetReader.h"
 #include "AliJetESDReader.h"
 #include "AliJetESDReaderHeader.h"
-#include "AliESD.h"
-//#include "AliEMCALGeometry.h"
+//#include "AliESD.h"
+#include "AliESDEvent.h"
 #include "AliJetDummyGeo.h"
 #include "AliJetUnitArray.h"
 #include "AliJetFillUnitArrayTracks.h"
@@ -77,59 +55,108 @@ ClassImp(AliJetFillUnitArrayTracks)
 
 //_____________________________________________________________________________
 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
-    : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
-      fNumUnits(0),
-      fEtaMinCal(0.),
-      fEtaMaxCal(0.),
-      fPhiMinCal(0.),
-      fPhiMaxCal(0.),
-      fHadCorr(0x0),
-      fHCorrection(0),
-      fNIn(0),
-      fOpt(0),
-      fDebug(0),
-      fReaderHeader(0x0),
-      fMomentumArray(0x0),
-      fUnitArray(0x0),
-      fTPCGrid(0x0),
-      fEMCalGrid(0x0),
-      fGeom(0x0),
-      fNphi(0),
-      fNeta(0),
-      fPhi2(0x0),
-      fEta2(0x0),
-      fPhi(0x0),
-      fEta(0x0),
-      fIndex(0x0),
-      fParams(0x0),
-      fGrid(0),
-      fPhiMin(0.),
-      fPhiMax(0.),
-      fEtaMin(0.),
-      fEtaMax(0.),
-      fEtaBinInTPCAcc(0),
-      fPhiBinInTPCAcc(0),
-      fEtaBinInEMCalAcc(0),
-      fPhiBinInEMCalAcc(0),
-      fNbinPhi(0)
+  : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
+    fNumUnits(0),
+    fEtaMinCal(0),
+    fEtaMaxCal(0),
+    fPhiMinCal(0),
+    fPhiMaxCal(0),
+    fHadCorr(0),
+    fHCorrection(0),
+    fNTracks(0),
+    fNTracksCut(0),
+    fOpt(0),
+    fDZ(0),
+    fDebug(0),
+    fReaderHeader(0x0),
+    fMomentumArray(0x0),
+    fUnitArray(0x0),
+    fRefArray(0x0),
+    fTPCGrid(0x0),
+    fEMCalGrid(0x0),
+    fGeom(0x0),
+    fESD(0x0),
+    fGrid0(0x0),
+    fGrid1(0x0),
+    fGrid2(0x0),
+    fGrid3(0x0),
+    fGrid4(0x0),
+    fNphi(0),
+    fNeta(0),
+    fPhi2(0x0),
+    fEta2(0x0),
+    fPhi(0x0),
+    fEta(0x0),
+    fIndex(0x0),
+    fParams(0x0),
+    fGrid(0),
+    fPhiMin(0),
+    fPhiMax(0),
+    fEtaMin(0),
+    fEtaMax(0),
+    fEtaBinInTPCAcc(0),
+    fPhiBinInTPCAcc(0),
+    fEtaBinInEMCalAcc(0),
+    fPhiBinInEMCalAcc(0),
+    fNbinPhi(0)
 {
   // constructor
 }
 
-//____________________________________________________________________________
-void AliJetFillUnitArrayTracks::SetEMCALGeometry()
+//_____________________________________________________________________________
+AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/)
+  : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
+    fNumUnits(0),
+    fEtaMinCal(0),
+    fEtaMaxCal(0),
+    fPhiMinCal(0),
+    fPhiMaxCal(0),
+    fHadCorr(0),
+    fHCorrection(0),
+    fNTracks(0),
+    fNTracksCut(0),
+    fOpt(0),
+    fDZ(0),
+    fDebug(0),
+    fReaderHeader(0x0),
+    fMomentumArray(0x0),
+    fUnitArray(0x0),
+    fRefArray(0x0),
+    fTPCGrid(0x0),
+    fEMCalGrid(0x0),
+    fGeom(0x0),
+    fESD(0x0),
+    fGrid0(0x0),
+    fGrid1(0x0),
+    fGrid2(0x0),
+    fGrid3(0x0),
+    fGrid4(0x0),
+    fNphi(0),
+    fNeta(0),
+    fPhi2(0x0),
+    fEta2(0x0),
+    fPhi(0x0),
+    fEta(0x0),
+    fIndex(0x0),
+    fParams(0x0),
+    fGrid(0),
+    fPhiMin(0),
+    fPhiMax(0),
+    fEtaMin(0),
+    fEtaMax(0),
+    fEtaBinInTPCAcc(0),
+    fPhiBinInTPCAcc(0),
+    fEtaBinInEMCalAcc(0),
+    fPhiBinInEMCalAcc(0),
+    fNbinPhi(0)
 {
-    // Set EMCAL geometry information
-    fGeom = AliJetDummyGeo::GetInstance();
-    if (fGeom == 0)
-       fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
-    if(fDebug>1) printf("\n EMCAL Geometry setted ! \n");
+  // constructor
 }
 
 //____________________________________________________________________________
 void AliJetFillUnitArrayTracks::InitParameters()
 {
-  fHCorrection    = 0;     // For hadron correction
+  //  fHCorrection    = 0;     // For hadron correction
   fHadCorr        = 0;     // For hadron correction
   fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
   fDebug = fReaderHeader->GetDebug();
@@ -139,13 +166,6 @@ void AliJetFillUnitArrayTracks::InitParameters()
   fPhiMinCal = fGeom->GetArm1PhiMin();
   fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
 
-  if(fDebug>30){
-    cout << "fEtaMinCal : " << fEtaMinCal << endl;
-    cout << "fEtaMaxCal : " << fEtaMaxCal << endl;
-    cout << "fPhiMinCal : " << fPhiMinCal << endl;
-    cout << "fPhiMaxCal : " << fPhiMaxCal << endl;
-  }
-
   fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, 
                        fPhiMax,fEtaMin,fEtaMax);
   fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc, 
@@ -173,349 +193,360 @@ AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
 }
 
 //_____________________________________________________________________________
-void AliJetFillUnitArrayTracks::Exec(Option_t* option)
+void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
 {
   //
   // Main method.
-  // Explain
+  //
 
   fDebug = fReaderHeader->GetDebug();
-  if(fDebug>1) printf("In AliJetFillUnitArrayTracks !");  
-  if(fDebug>3) printf("\nfGeom->GetEntries() = %d\n", fGeom->GetNCells());
-  // Set EMCal Geometry
-  SetEMCALGeometry();
+
   // Set parameters
   InitParameters();
 
-  TClonesArray *lvArray = fMomentumArray;     // Correct checked !
-  Int_t nInT =  lvArray->GetEntries();        // Correct checked !
-  Float_t ptMin = fReaderHeader->GetPtCut();  // Correct checked !
-  
-  // sflag -> not yet implemented !!!!
-  
-  if(fDebug>3) cout << "nInT : " << nInT << endl;
-
-  if (nInT == 0) return;
-  
-  // local arrays for input
-  Float_t* enT  = new Float_t[nInT];
-  Float_t* ptT  = new Float_t[nInT];
-  Float_t* etaT = new Float_t[nInT];
-  Float_t* phiT = new Float_t[nInT];
-  Float_t* thetaT = new Float_t[nInT];
-  
-  Int_t trackInEmcalAcc = 0;
-  Int_t trackInTpcAcc = 0;
-  Int_t trackInTpcAccOnly = 0;
-
-  Int_t nElements = fTPCGrid->GetNEntries();
-  Int_t nElements2 = fEMCalGrid->GetNEntries();
-  fGrid = fTPCGrid->GetGridType();
+  // get number of tracks in event (for the loop)
+  Int_t goodTrack = 0;
+  Int_t nt = 0;
+  Float_t pt, eta,phi;
+  TVector3 p3;
+  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();
 
-  if(fDebug>3){
-    cout << "nElements : " << nElements << endl;
-    cout << "nElements2 : " << nElements2 << endl;
-    cout << "fNumUnits : " << fNumUnits << endl;
-    cout << "sum : " << fNumUnits+nElements << endl;
-  }
+  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;
 
-  // Set energy exactly to zero
-  if(fGrid==0)
-    for(Int_t k=0; k<nElements; k++)
-      fUnitArray[k].SetUnitEnergy(0.); 
-  
-  if(fGrid==1)
-    for(Int_t k=0; k<fNumUnits+nElements; k++)
-      fUnitArray[k].SetUnitEnergy(0.);  
-
-  // load input vectors
-  for (Int_t i = 0; i < nInT; i++)  
-    {
-      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
-      enT[i]  = lv->Energy();
-      ptT[i]  = lv->Pt();
-      etaT[i] = lv->Eta();
-      phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
-      thetaT[i] = 2.0*TMath::ATan(TMath::Exp(-etaT[i]));
-      
-      if(fDebug>20){    
-       cout << "enT[" << i << "] : " <<  enT[i] << endl;
-       cout << "ptT[" << i << "] : " <<  ptT[i] << endl;
-       cout << "etaT[" << i << "] : " <<  etaT[i] << endl;
-       cout << "phiT[" << i << "] : " <<  phiT[i] << endl;
-       cout << "thetaT[" << i << "] : " <<  thetaT[i] << endl;
-       cout << "fEtaMinCal : " << fEtaMinCal << ", fEtaMaxCal : " << fEtaMaxCal << endl;
-       cout << "fPhiMinCal : " << fPhiMinCal << ", fPhiMaxCal : " << fPhiMaxCal << endl;
-       cout << "fEtaMin : " << fEtaMin << ", fEtaMax : " << fEtaMax << endl;
-       cout << "fPhiMin : " << fPhiMin << ", fPhiMax : " << fPhiMax << endl;
-      }
-      
-      if(fGrid==0)
-       {
-         // For the moment, only TPC filled from its grid in its total acceptance
-         if(fDebug>2)
-           cout << "In total TPC acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
-           
-           trackInTpcAccOnly += 1;
-           
-           Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
-           
-           Float_t unitEnergy = 0.;
-           unitEnergy = fUnitArray[idTPC].GetUnitEnergy();
-           
-           if(unitEnergy > 0. && fDebug >10){
-             cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-             cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-             cout << "TPC id : " << idTPC << endl;
-             cout << "etaT[" << i << "] : " << etaT[i] << endl;
-             cout << "phiT[" << i << "] : " << phiT[i] << endl;
-             cout << "unitEnergy in TPC acceptance : " << unitEnergy << endl;
-             cout << "fUnitArray[idTPC].GetUnitEnergy(): " << 
-               fUnitArray[idTPC].GetUnitEnergy() << endl;
-             cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-             cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-           }
+  fGrid = fTPCGrid->GetGridType();
 
-           // Fill energy in TPC acceptance
-           fUnitArray[idTPC].SetUnitEnergy(unitEnergy + ptT[i]);
 
-           if(fDebug > 10){
-             cout << "ptT[" << i << "] : " << ptT[i] << endl;
-             cout << "unitEnergy in TPC acceptance after : " << 
-               fUnitArray[idTPC].GetUnitEnergy() << endl;
-           }
-           
-           // Pt cut flag
-           if(fUnitArray[idTPC].GetUnitEnergy()<ptMin)
-             fUnitArray[idTPC].SetUnitCutFlag(kPtSmaller);
-           else fUnitArray[idTPC].SetUnitCutFlag(kPtHigher);
-           // Detector flag
-           fUnitArray[idTPC].SetUnitDetectorFlag(kTpc);
+  //loop over tracks
+  for (Int_t it = 0; it < nt; it++) {
+    AliESDtrack *track = fESD->GetTrack(it);
+    UInt_t status = track->GetStatus();
+    
+    Double_t mom[3];
+    track->GetPxPyPz(mom);
+    p3.SetXYZ(mom[0],mom[1],mom[2]);
+    pt = p3.Pt();
+    Float_t mass = 0.;
+    mass = track->GetMass();
+    
+    if (((status & AliESDtrack::kITSrefit) == 0) ||
+        ((status & AliESDtrack::kTPCrefit) == 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 -> not yet implemented !!!!
+    
+    if(fGrid==0)
+      {
+       // Only TPC filled from its grid in its total acceptance
+       
+       Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
+       Bool_t ok = kFALSE;
+
+       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
+       uArray->SetUnitTrackID(it);
+
+       Float_t unitEnergy = 0.;
+       unitEnergy = uArray->GetUnitEnergy();
+       if(unitEnergy==0.){
+         nTracksTpcOnly++;
+         ok = kTRUE;
+       }
+       // Fill energy in TPC acceptance
+       uArray->SetUnitEnergy(unitEnergy + pt);
+       uArray->SetUnitPxPyPz(mom);
+       uArray->SetUnitMass(mass);
+
+       // Pt cut flag
+       if(uArray->GetUnitEnergy()<ptMin){
+         uArray->SetUnitCutFlag(kPtSmaller);
+       }
+       else {
+         uArray->SetUnitCutFlag(kPtHigher);
+         if(ok) nTracksTpcOnlyCut++;
        }
 
-      if(fGrid==1)
-       {
-         // Fill track information in EMCAL acceptance
-         if((etaT[i] >= fEtaMin && etaT[i] <= fEtaMax) &&
-            (phiT[i] >= fPhiMin && phiT[i] <= fPhiMax))// &&
-           //   GetCutFlag(i) == 1)
-           //   ptT[i] > ptMin)
-           {
-             trackInEmcalAcc += 1;
-             
-             if(fDebug>20){ 
-               cout << "before : " << endl;
-               cout << "etaT[i] : " << etaT[i] << endl;
-               cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
-             }
-             
-             // This function should be modified soon
-//           Int_t towerID = fGeom->TowerIndexFromEtaPhi(etaT[i],180.0/TMath::Pi()*phiT[i]); // Obsolete
-             Int_t towerID = fGeom->TowerIndexFromEtaPhi2(etaT[i],180.0/TMath::Pi()*phiT[i]); // Mine modified
-//           Int_t towerID = fEMCalGrid->GetIndexFromEtaPhi(phiT[i],etaT[i]);  // Using an EMCal grid -> calculated
-//           Int_t towerID = fEMCalGrid->GetIndex(phiT[i],etaT[i]);  // Using an EMCal grid -> tabulated (faster)
-
-             Float_t unitEnergy = fUnitArray[towerID].GetUnitEnergy(); 
-             
-             if(fDebug>20) { 
-               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-               cout << "after : " << endl;
-               cout << "towerID : " << towerID << endl;
-               cout << "etaT[i] : " << etaT[i] << endl;
-               cout << "phiT[i](rad) : " << phiT[i] << endl;
-               cout << "phiT[i] : " << phiT[i]*180/TMath::Pi() << endl;
-               cout << "unitEnergy in emcal acceptance : " << unitEnergy << endl;
-               cout << "fHadCorr : " << fHadCorr << endl;
-               cout << "fHCorrection : " << fHCorrection << endl;
-               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-               cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-             }
-
-             //OLD WAY:   //Do Hadron Correction
-             if (fHCorrection != 0) 
-               { 
-                 Float_t   hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
-                 unitEnergy -= hCEnergy*TMath::Sin(thetaT[i]);
-                 if(fDebug>20){
-                   cout << "Inside loop for hadron correction ==========================" << endl;
-                   cout << "enT[" << i << "] : " << enT[i] << endl;
-                   cout << "ptT[" << i << "] : " << ptT[i] << endl;
-                   cout << "etaT[" << i << "] : " << etaT[i] << endl;
-                   cout << "Energy correction : " << hCEnergy << endl;
-                   cout << "unitEnergy : " << unitEnergy << endl;
-                   cout << "fUnitArray[towerID].GetUnitEnergy() : " << 
-                     fUnitArray[towerID].GetUnitEnergy() << endl;
-                 }           
-               } //end Hadron Correction loop
-             
-             fUnitArray[towerID].SetUnitEnergy(unitEnergy + ptT[i]);
-
-             // Put a pt cut flag
-             if(fUnitArray[towerID].GetUnitEnergy()<ptMin)
-               fUnitArray[towerID].SetUnitCutFlag(kPtSmaller);
-             else fUnitArray[towerID].SetUnitCutFlag(kPtHigher);
-             // Detector flag
-             fUnitArray[towerID].SetUnitDetectorFlag(kTpc);
-             
-             if(fDebug>10){
-               cout << "After pT filled ===============================================" << endl;
-               cout << "PtCut : " << ptMin << endl;
-               cout << "ptT[" << i << "] : " << ptT[i] << endl;
-               cout << "unitEnergy : " << unitEnergy << endl;
-               cout << "fUnitArray[towerID].GetUnitEnergy() : " << fUnitArray[towerID].GetUnitEnergy() << endl;
-               cout << "fUnitArray[towerID].GetUnitCutFlag() : " << fUnitArray[towerID].GetUnitCutFlag() << endl;
-               cout << "fUnitArray[towerID].GetUnitEta() : " << fUnitArray[towerID].GetUnitEta() << endl;
-               cout << "fUnitArray[towerID].GetUnitPhi() : " << fUnitArray[towerID].GetUnitPhi() << endl;
-             }
-           } // end loop on EMCal acceptance cut + tracks quality
-         else{ 
-           // Outside EMCal acceptance
-           if(fDebug>2)
-             cout << "Outside EMCal acceptance +++++++++++++++++++++++++++++++++++++++++++" << endl;
-           
-           trackInTpcAcc += 1;
-           
-           //  Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
-           Int_t idTPC = fTPCGrid->GetIndex(phiT[i],etaT[i]);
-           
-           Float_t unitEnergy2 = 0.;
-           unitEnergy2 = fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
-           
-           if(unitEnergy2 > 0. && fDebug >10){
-             cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-             cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-             cout << "TPC id : " << idTPC << endl;
-             cout << "Total id : " << fNumUnits-1+idTPC << endl;
-             cout << "etaT[" << i << "] : " << etaT[i] << endl;
-             cout << "phiT[" << i << "] : " << phiT[i] << endl;
-             cout << "unitEnergy outside emcal acceptance : " << unitEnergy2 << endl;
-             cout << "fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy(): " << 
-               fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy() << endl;
-             cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-             cout << "§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§§" << endl;
-           }
+       // Detector flag
+       if(unitEnergy>0) {
+         uArray->SetUnitDetectorFlag(kAll);
+       }
+       if(uArray->GetUnitEnergy()>0){
+         fRefArray->Add(uArray);
+       }
 
-           // Fill energy outside emcal acceptance
-           fUnitArray[fNumUnits-1+idTPC].SetUnitEnergy(unitEnergy2 + ptT[i]);
+       sflag[goodTrack]=0;
+       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+       cflag[goodTrack]=0;
+       if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
+       goodTrack++;
 
-           // Pt cut flag
-           if(fUnitArray[fNumUnits-1+idTPC].GetUnitEnergy()<ptMin)
-             fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtSmaller);
-           else fUnitArray[fNumUnits-1+idTPC].SetUnitCutFlag(kPtHigher);
-           // Detector flag
-           fUnitArray[fNumUnits-1+idTPC].SetUnitDetectorFlag(kTpc);
-         }
-       } // end fGrid==1
-    } // end loop on entries (tpc tracks)
-      
-  if(fDebug>3){
-    printf("Number of tracks in EMCAL acceptance: %d\n", trackInEmcalAcc);
-    printf("Number of tracks in TPC acceptance: %d\n", trackInTpcAcc);
-  }
-  
-  if(fGrid==0) fNIn = trackInTpcAccOnly;
-  if(fGrid==1) fNIn = trackInEmcalAcc+trackInTpcAcc;
-  
-  fOpt = fReaderHeader->GetDetector();
-  
-  if(fDebug>1) printf("fOpt in FillUnitArrayFromTPCTracks : %d\n", fOpt);
-  
-      if(fOpt==1 || option=="tpc") // if only TPC
-       { //Set all unit flags, Eta, Phi
-
-         if(fGrid==0)
-           {
-             if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
-             for(Int_t i=0; i<nElements; i++)
-               {
-                 Float_t eta = 0.;
-                 Float_t phi = 0.;
-                 if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
-                 fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
-                 fUnitArray[i].SetUnitEntries(nElements);
-                 fUnitArray[i].SetUnitID(i);
-                 fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta);
-                 fUnitArray[i].SetUnitEta(eta);
-                 fUnitArray[i].SetUnitPhi(phi);
-                 fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
-                 fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
+      }
+    
+    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();
+
+               if(phi >= phimin0 && phi <= phimax0){
+                 Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
+                 AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
+                 uArray0->SetUnitTrackID(it);
+                 Float_t uEnergy0 = uArray0->GetUnitEnergy();
+                 Bool_t ok0 = kFALSE;
+                 if(uEnergy0==0.){
+                   nTracksEmcalDZ++;
+                   ok0 = kTRUE;
+                 }
+                 uArray0->SetUnitEnergy(uEnergy0+pt);
+                 if(uArray0->GetUnitEnergy()<ptMin)
+                   uArray0->SetUnitCutFlag(kPtSmaller);
+                 else {
+                   uArray0->SetUnitCutFlag(kPtHigher);
+                   if(ok0) nTracksEmcalDZCut++;
+                 }
+                 if(uArray0->GetUnitEnergy()>0)
+                   fRefArray->Add(uArray0);
                }
-           }
-
-         if(fGrid==1)
-           {
-             if(fDebug>1) printf("In if(fOpt==1 || option==tpc)\n");
-             for(Int_t i=0; i<fNumUnits+nElements; i++)
-               {
-                 Float_t eta = 0.;
-                 Float_t phi = 0.;
-                 if (fDebug>10) Info("FillUnitArray","Setting all units outside jets \n");
-                 //Set all units to be outside a jet initially
-                 fUnitArray[i].SetUnitFlag(kOutJet); // returns 0, 1, 2...
-                 fUnitArray[i].SetUnitEntries(fNumUnits+nElements);
-                 fUnitArray[i].SetUnitID(i);
-
-                 if(fUnitArray[i].GetUnitID()<fNumUnits)
-                   {
-                     // fGeom->EtaPhiFromIndex2(fUnitArray[i].GetUnitID(), eta, phi); // My function in HEADPythia63
-                     fGeom->EtaPhiFromIndex(fUnitArray[i].GetUnitID(), eta, phi); // From EMCal geometry 
-                     phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
-                     // fEMCalGrid->GetEtaPhiFromIndex2(i,phi,eta); // My function from Grid
-                     // fEMCalGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID(),phi,eta); // My function from Grid
-                     fUnitArray[i].SetUnitEta(eta);
-                     //fUnitArray[i].SetUnitPhi(phi*TMath::Pi()/180.0);
-                     fUnitArray[i].SetUnitPhi(phi);
-                     fUnitArray[i].SetUnitDeta(fEMCalGrid->GetDeta());
-                     fUnitArray[i].SetUnitDphi(fEMCalGrid->GetDphi());
-                   } 
+               if(phi >= phimin1 && phi <= phimax1){
+                 Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
+                 AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
+                 uArray1->SetUnitTrackID(it);
+                 Float_t uEnergy1 = uArray1->GetUnitEnergy();
+                 Bool_t ok1 = kFALSE;
+                 if(uEnergy1==0.){
+                   nTracksEmcalDZ++;
+                   ok1 = kTRUE;
+                 }
+                 uArray1->SetUnitEnergy(uEnergy1+pt);
+                 if(uArray1->GetUnitEnergy()<ptMin)
+                   uArray1->SetUnitCutFlag(kPtSmaller);
                  else {
-                   fTPCGrid->GetEtaPhiFromIndex2(fUnitArray[i].GetUnitID()+1-fNumUnits,phi,eta);
-                   fUnitArray[i].SetUnitEta(eta);
-                   fUnitArray[i].SetUnitPhi(phi);
-                   fUnitArray[i].SetUnitDeta(fTPCGrid->GetDeta());
-                   fUnitArray[i].SetUnitDphi(fTPCGrid->GetDphi());
-                   
-                   if(fDebug>10)
-                     {
-                       if(fUnitArray[i].GetUnitEnergy()!=0.){
-                         cout << "(fUnitArray[" << i << "].GetUnitID()+1-fNumUnits : " << 
-                           fUnitArray[i].GetUnitID()+1-fNumUnits << endl;
-                         cout << "(fUnitArray[" << i << "].GetUnitEnergy() : " << 
-                           fUnitArray[i].GetUnitEnergy() << endl;
-                         cout << "(fUnitArray[" << i << "].GetUnitEta() : " << 
-                           fUnitArray[i].GetUnitEta() << endl;
-                         cout << "(fUnitArray[" << i << "].GetUnitPhi() : " << 
-                           fUnitArray[i].GetUnitPhi() << endl;
-                       }
-                     } 
+                   uArray1->SetUnitCutFlag(kPtHigher);
+                   if(ok1) nTracksEmcalDZCut++;
                  }
-                 fUnitArray[i].SetUnitClusterID(0);
-               }//end loop over all units in array (same as all towers in EMCAL)
-           }
-
-         if(fDebug>20) 
-           {
-             for(Int_t i=0; i<fNumUnits+nElements; i++)
-               {
-                 if(fUnitArray[i].GetUnitEnergy()!=0.){
-                   cout << "######################################################### " << endl;
-                   cout << "Final UnitArray filled with energy != 0" << i << endl;
-                   cout << "Pointeur UnitArray : " << fUnitArray << " ID : " << 
-                     fUnitArray[i].GetUnitID() << " Energy : " << fUnitArray[i].GetUnitEnergy() << 
-                     " Eta : " << fUnitArray[i].GetUnitEta() << " Phi : " << fUnitArray[i].GetUnitPhi() << endl;
+                 if(uArray1->GetUnitEnergy()>0)
+                   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);
+                 uArray2->SetUnitTrackID(it);
+                 Float_t uEnergy2 = uArray2->GetUnitEnergy();
+                 Bool_t ok2 = kFALSE;
+                 if(uEnergy2==0.){
+                   nTracksEmcalDZ++;
+                   ok2 = kTRUE;
+                 }
+                 uArray2->SetUnitEnergy(uEnergy2+pt);
+                 if(uArray2->GetUnitEnergy()<ptMin)
+                   uArray2->SetUnitCutFlag(kPtSmaller);
+                 else {
+                   uArray2->SetUnitCutFlag(kPtHigher);
+                   if(ok2) nTracksEmcalDZCut++;
+                 }
+                 if(uArray2->GetUnitEnergy()>0)
+                   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);
+                 uArray3->SetUnitTrackID(it);
+                 Float_t uEnergy3 = uArray3->GetUnitEnergy();
+                 Bool_t ok3 = kFALSE;
+                 if(uEnergy3==0.){
+                   nTracksEmcalDZ++;
+                   ok3 = kTRUE;
+                 }
+                 uArray3->SetUnitEnergy(uEnergy3+pt);
+                 if(uArray3->GetUnitEnergy()<ptMin)
+                   uArray3->SetUnitCutFlag(kPtSmaller);
+                 else {
+                   uArray3->SetUnitCutFlag(kPtHigher);
+                   if(ok3) nTracksEmcalDZCut++;
+                 }
+                 if(uArray3->GetUnitEnergy()>0)
+                   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);
+                 uArray4->SetUnitTrackID(it);
+                 Float_t uEnergy4 = uArray4->GetUnitEnergy();
+                 Bool_t ok4 = kFALSE;
+                 if(uEnergy4==0.){
+                   nTracksEmcalDZ++;
+                   ok4 = kTRUE;
+                 }
+                 uArray4->SetUnitEnergy(uEnergy4+pt);
+                 if(uArray4->GetUnitEnergy()<ptMin)
+                   uArray4->SetUnitCutFlag(kPtSmaller);
+                 else {
+                   uArray4->SetUnitCutFlag(kPtHigher);
+                   if(ok4) nTracksEmcalDZCut++;
                  }
+                 if(uArray4->GetUnitEnergy()>0)
+                   fRefArray->Add(uArray4);
                }
+             } // end fDZ
+           
+           Int_t towerID = 0;
+           fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
+
+           if(towerID==-1) continue;
+
+           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
+           uArray->SetUnitTrackID(it);
+           Float_t unitEnergy = uArray->GetUnitEnergy(); 
+           Bool_t ok = kFALSE;
+           if(unitEnergy==0.){
+             nTracksEmcal++;
+             ok=kTRUE;
+           }
+
+           // Do Hadron Correction
+           // Parametrization to be added
+           if (fHCorrection != 0) 
+             { 
+               //            Float_t   hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
+               Float_t   hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta);
+               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++;
            }
+           // Detector flag
+           if(unitEnergy > 0)
+             uArray->SetUnitDetectorFlag(kAll);
+           
+           if(uArray->GetUnitEnergy()>0)
+             fRefArray->Add(uArray);
+               
+           sflag[goodTrack]=0;
+           if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+           cflag[goodTrack]=0;
+           if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
+           goodTrack++;
+           
+         } // 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);
+         uArray->SetUnitTrackID(it);
+
+         Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
+         Bool_t ok2 = kFALSE;
+         if(unitEnergy2==0.){
+           nTracksTpc++;
+           ok2=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(ok2) nTracksTpcCut++;
+         }
+         // Detector flag
+         if(unitEnergy2 > 0)
+           uArray->SetUnitDetectorFlag(kTpc);
+         if(uArray->GetUnitEnergy()>0)
+           fRefArray->Add(uArray);
+
+         sflag[goodTrack]=0;
+         if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+         cflag[goodTrack]=0;
+         if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
+         goodTrack++;
+         
+       }
+      } // end fGrid==1
+  } // end loop on entries (tpc tracks)
 
-       } // end  if(fOpt==1 || option=="tpc")
+//   // set the signal flags
+//   fSignalFlag.Set(goodTrack,sflag);
+//   fCutFlag.Set(goodTrack,cflag);
 
-      delete[] enT;
-      delete[] ptT;
-      delete[] etaT;
-      delete[] phiT;
-      delete[] thetaT;
+//   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;
+    }
+  }  
+
 }
 
 //__________________________________________________________
index 08e6356..42fc6ea 100644 (file)
 #include <TMatrixD.h>
 #include <TArrayD.h>
 
-class AliEMCALGeometry;
-class AliJetDummyGeo;
 class AliJetHadronCorrection;
 class AliJetReader;
 class AliJetESDReader;
 class TClonesArray;
-class AliJetUnitArray;
+class TRefArray;
+//class AliJetUnitArray;
 //class AliJetReaderHeader;
-class AliJetReader;
 class AliJetGrid;
+class AliJetDummyGeo;
+class AliESD;
+class AliESDEvent;
 
 class AliJetFillUnitArrayTracks : public TTask
 {
  public: 
   AliJetFillUnitArrayTracks();
+  AliJetFillUnitArrayTracks(AliESD *fESD);
+  AliJetFillUnitArrayTracks(AliESDEvent *fESD);
   virtual ~AliJetFillUnitArrayTracks();
   
   // Setter
   void SetReaderHeader(AliJetReaderHeader *readerHeader) {fReaderHeader = readerHeader;}
+  void SetGeom(AliJetDummyGeo *geom) {fGeom = geom;}
   void SetMomentumArray(TClonesArray *momentumArray) {fMomentumArray = momentumArray;}
-  void SetUnitArray(AliJetUnitArray *unitArray) {fUnitArray = unitArray;}
+  void SetUnitArray(TClonesArray *unitArray) {fUnitArray = unitArray;}
+  void SetRefArray(TRefArray *refArray) {fRefArray = refArray;}
   void SetHadCorrection(Int_t flag = 1) {fHCorrection = flag;}
   void SetHadCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;}
   void SetTPCGrid(AliJetGrid *grid) {fTPCGrid = grid;}
   void SetEMCalGrid(AliJetGrid *grid) {fEMCalGrid = grid;}
   void SetGrid(Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax);
+  //  void SetESD(AliESD *esd) {fESD = esd;}
+  void SetESD(AliESDEvent *esd) {fESD = esd;}
+  void SetGrid0(AliJetGrid *grid0){fGrid0 = grid0;}
+  void SetGrid1(AliJetGrid *grid1){fGrid1 = grid1;}
+  void SetGrid2(AliJetGrid *grid2){fGrid2 = grid2;}
+  void SetGrid3(AliJetGrid *grid3){fGrid3 = grid3;}
+  void SetGrid4(AliJetGrid *grid4){fGrid4 = grid4;}
 
   // Getter
-  AliJetUnitArray* GetUnitArray() {return fUnitArray;}
-  //  Int_t GetIndexFromEtaPhi(Double_t eta,Double_t phi) const;
-  void  GetEtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi);
-  Int_t GetNeta() {return fNeta;}
-  Int_t GetNphi() {return fNphi;}
-
-  void Exec(Option_t*);
+  TClonesArray* GetUnitArray() {return fUnitArray;}
+  TRefArray*    GetRefArray() {return fRefArray;}
+  //  Int_t         GetIndexFromEtaPhi(Double_t eta,Double_t phi) const;
+  void          GetEtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi);
+  Int_t         GetNeta() {return fNeta;}
+  Int_t         GetNphi() {return fNphi;}
+  Int_t         GetHadCorrection() {return fHCorrection;}
+  Int_t         GetMult() {return fNTracks;}
+  Int_t         GetMultCut() {return fNTracksCut;}
+  void          Exec(Option_t*);
 
  protected:
   Int_t   fNumUnits;      // Number of units in the unit object array (same as num towers in EMCAL)
@@ -60,19 +75,27 @@ class AliJetFillUnitArrayTracks : public TTask
   Float_t fPhiMinCal;     // Define EMCal acceptance in Phi
   Float_t fPhiMaxCal;     // Define EMCal acceptance in Phi
   AliJetHadronCorrectionv1   *fHadCorr;         // Pointer to Hadron Correction Object
-  Int_t                       fHCorrection;     // Hadron correction flag
-  Int_t                       fNIn;             // Number of Array filled in UnitArray
+  Int_t                       fHCorrection;     //  Hadron correction flag
+  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
-  AliJetUnitArray             *fUnitArray;      // UnitArray
+  TClonesArray                *fUnitArray;      // UnitArray
+  TRefArray                   *fRefArray;       // UnitArray
   AliJetGrid                  *fTPCGrid;        // Define filled grid
   AliJetGrid                  *fEMCalGrid;      // Define filled grid
-
-  // geometry info
-  AliJetDummyGeo              *fGeom;           //!
+  AliJetDummyGeo              *fGeom;           // Define EMCal geometry
+  //  AliESD                     *fESD;            // ESD
+  AliESDEvent                 *fESD;            // ESD
+  AliJetGrid                  *fGrid0;
+  AliJetGrid                  *fGrid1;
+  AliJetGrid                  *fGrid2;
+  AliJetGrid                  *fGrid3;
+  AliJetGrid                  *fGrid4;
 
   Int_t     fNphi;                    // number of points in the grid:   phi
   Int_t     fNeta;                    //               "                 eta
@@ -95,11 +118,9 @@ class AliJetFillUnitArrayTracks : public TTask
   Int_t     fNbinPhi;
 
  private:
-  void SetEMCALGeometry();
+  //  void SetEMCALGeometry();
   void InitParameters();
 
-  //  Int_t fEvent;
-  
   ClassDef(AliJetFillUnitArrayTracks,1) // Fill Unit Array with tpc and/or emcal information
 };
 
index 58d19fb..02416cc 100644 (file)
@@ -253,7 +253,5 @@ void AliJetFinder::AddJet(AliAODJet p)
 void AliJetFinder::ConnectAOD(AliAODEvent* aod)
 {
 // Connect to the AOD
-    printf("Connect AOD \n");
     fAODjets = aod->GetJets();
-    printf("Connect AOD %p \n", fAODjets);
 }
index 61503d0..5ac309a 100644 (file)
@@ -26,7 +26,6 @@
 // --- AliRoot header files ---
 //#include "../EMCAL/AliJetGeometry.h"
 #include "AliJetDummyGeo.h"
-//#include "AliEMCALGeometry.h"
 #include "AliJetHadronCorrectionv1.h"
 
 ClassImp(AliJetHadronCorrectionv1)
index b8da7fc..36dc663 100644 (file)
@@ -11,7 +11,6 @@
 #define HCPARAMETERS    6
 #define HCPARAMETERSETS 2
 
-class AliEMCALGeometry;
 class AliJetDummyGeo;
 
 
index 6007862..7ebb8ac 100644 (file)
@@ -25,6 +25,7 @@
 #include <TVector3.h>
 #include <TLorentzVector.h>
 #include <TSystem.h>
+#include <TChain.h>
 // From AliRoot ...
 #include "AliJetMCReader.h"
 #include "AliJetMCReaderHeader.h"
index 2f6dd92..db62f0f 100755 (executable)
 
 // root
 #include <TClonesArray.h>
+#include <TRefArray.h>
+#include "TTask.h"
 //AliRoot
 #include "AliJetReader.h"
 #include "AliJetReaderHeader.h"
+#include "AliESDEvent.h"
+#include "AliHeader.h"
+#include "AliJetFillUnitArrayTracks.h" 
+#include "AliJetFillUnitArrayEMCalDigits.h"
 #include "AliJetUnitArray.h"
 #include "AliJetHadronCorrectionv1.h"
 
@@ -33,15 +39,23 @@ ClassImp(AliJetReader)
 ////////////////////////////////////////////////////////////////////////
 
 AliJetReader::AliJetReader():
+  // Constructor
+  fChain(0), 
   fMomentumArray(new TClonesArray("TLorentzVector",2000)),
   fArrayMC(0),
   fFillUnitArray(new TTask("fillUnitArray","Fill unit array jet finder")),
+  fESD(0),
   fReaderHeader(0),
   fSignalFlag(0),
   fCutFlag(0),
-  fUnitArray(new AliJetUnitArray[60000]),     
-  fUnitArrayNoCuts(new AliJetUnitArray[60000]),
-  fArrayInitialised(0)
+  fUnitArray(new TClonesArray("AliJetUnitArray",60000)),
+  fRefArray(new TRefArray()),
+  fUnitArrayNoCuts(new TClonesArray("AliJetUnitArray",60000)),
+  fArrayInitialised(0),
+  fFillUAFromTracks(new AliJetFillUnitArrayTracks()), 
+  fFillUAFromEMCalDigits(new AliJetFillUnitArrayEMCalDigits())
+  //  fHadronCorrector(0),
+  //  fHCorrection(0)
 {
   // Default constructor
   fSignalFlag = TArrayI();
index 3b63c32..a104c34 100755 (executable)
 #include <TObject.h>
 #include <TChain.h>
 #include <TArrayI.h>
-#ifndef ROOT_TTask
-#include "TTask.h"
-#endif
 
 class TTree;
 class TTask;
+class TChain;
 class TClonesArray;
 class TRefArray;
 class AliJetReaderHeader;
+class AliESDEvent;
+class AliHeader;
 class AliJetUnitArray;
 class AliJetHadronCorrectionv1;
 class AliJet;
+class AliJetFillUnitArrayTracks;
+class AliJetFillUnitArrayEMCalDigits;
 
 class AliJetReader : public TObject 
 {
@@ -32,15 +34,18 @@ class AliJetReader : public TObject
   virtual ~AliJetReader();
 
   // Getters
-  virtual TClonesArray*        GetMomentumArray()   const {return fMomentumArray;}
-  virtual TRefArray*           GetReferences()      const {return 0;}   
-  virtual AliJetUnitArray*     GetUnitArray()       const {return fUnitArray;}  
-  virtual AliJetUnitArray*     GetUnitArrayNoCuts() const {return fUnitArrayNoCuts;}
-  
-  virtual AliJetReaderHeader*  GetReaderHeader()    const {return fReaderHeader;}
-  virtual Int_t GetSignalFlag(Int_t i)              const {return fSignalFlag[i];}
-  virtual Int_t GetCutFlag(Int_t i)                 const {return fCutFlag[i];}
-  virtual Int_t GetArrayInitialised()               const {return fArrayInitialised;}
+  virtual TClonesArray*        GetMomentumArray()    const {return fMomentumArray;}
+  virtual TRefArray*           GetReferences()       const {return 0;}   
+  virtual TClonesArray        *GetUnitArray()        const {return fUnitArray;}  
+  virtual TRefArray           *GetRefArray()         const {return fRefArray;}
+  virtual TClonesArray        *GetUnitArrayNoCuts()  const {return fUnitArrayNoCuts;} 
+
+  virtual AliJetReaderHeader* GetReaderHeader()      const {return fReaderHeader;}
+  virtual Int_t               GetSignalFlag(Int_t i) const {return fSignalFlag[i];}
+  virtual Int_t               GetCutFlag(Int_t i)    const {return fCutFlag[i];}
+  virtual Int_t               GetArrayInitialised()  const {return fArrayInitialised;}
+  virtual Int_t               GetNumCandidate()      const {return fNumCandidate;} 
+  virtual Int_t               GetNumCandidateCut()   const {return fNumCandidateCut;}
   
   // Setters
   virtual Bool_t FillMomentumArray(Int_t) {return kTRUE;}
@@ -49,9 +54,20 @@ class AliJetReader : public TObject
   virtual void   FillUnitArrayFromEMCALDigits(Int_t) {}   // temporarily not used
   virtual void   FillUnitArrayFromEMCALClusters(Int_t) {} // temporarily not used
   virtual void   InitUnitArray() {}
+  virtual void   InitParameters() {}
+  virtual void   CreateTasks() {}
+  //  virtual void   ExecTasks(Int_t) {}
+  virtual Bool_t   ExecTasks(Int_t) {return kTRUE;}
+  /*   // Correction of hadronic energy 
+       virtual void SetHadronCorrector(AliEMCALHadronCorrectionv1* corr) {fHadronCorrector = corr;} 
+       virtual void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;} */
   virtual void   SetReaderHeader(AliJetReaderHeader* header) 
       {fReaderHeader = header;}
-         
+  virtual void   SetESD(AliESDEvent* esd) { fESD = esd;}
+  //  virtual Int_t  SetNumCandidate(Int_t cand) {fNumCandidate = cand;} 
+  //  virtual Int_t  SetNumCandidateCut(Int_t candcut) {fNumCandidateCut = candcut;}
+  
+
   // Others
   virtual void   OpenInputFiles() {}
   virtual void   SetInputEvent(TObject* /*esd*/, TObject* /*aod*/, TObject* /*mc*/) {;}
@@ -63,16 +79,23 @@ class AliJetReader : public TObject
  protected:
   AliJetReader(const AliJetReader& rJetReader);
   AliJetReader& operator = (const AliJetReader& rhsr);
+  TChain                  *fChain;            // chain for reconstructed tracks
   TClonesArray            *fMomentumArray;    // array of particle momenta
   TClonesArray            *fArrayMC;          //! array of mc particles
   TTask                   *fFillUnitArray;    //! task list for filling the UnitArray
+  AliESDEvent             *fESD;              // pointer to esd
   AliJetReaderHeader      *fReaderHeader;     // pointer to header
   TArrayI                  fSignalFlag;       // to flag if a particle comes from pythia or
                                               // from the underlying event
   TArrayI                  fCutFlag;          // to flag if a particle passed the pt cut or not
-  AliJetUnitArray         *fUnitArray;        //! array of digit position and energy 
-  AliJetUnitArray         *fUnitArrayNoCuts;  //! array of digit position and energy 
+  TClonesArray            *fUnitArray;        // array of digit position and energy 
+  TRefArray               *fRefArray;         // 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;
+  AliJetFillUnitArrayEMCalDigits   *fFillUAFromEMCalDigits;
+  Int_t                    fNumCandidate;
+  Int_t                    fNumCandidateCut;
   ClassDef(AliJetReader,1)
 };
  
index f172d84..7642a8d 100755 (executable)
@@ -19,6 +19,8 @@
 // Author: jgcn@mda.cinvestav.mx
 //-------------------------------------------------------------------------
 
+#include <TMath.h>
+
 #include "AliJetReaderHeader.h"
 
 ClassImp(AliJetReaderHeader)
@@ -30,9 +32,12 @@ AliJetReaderHeader::AliJetReaderHeader():
  fFirst(0),
  fLast(-1),
  fOption(0),
+ fDZ(0),
  fSignalPerBg(0),
  fFiducialEtaMin(-0.9),
  fFiducialEtaMax(0.9),
+ fFiducialPhiMin(0.),
+ fFiducialPhiMax(2*TMath::Pi()),
  fPtCut(2.0),
  fComment("No comment"),
  fDir(""),
@@ -49,9 +54,12 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name):
  fFirst(0),
  fLast(-1),
  fOption(0),
+ fDebug(0),
  fSignalPerBg(0),
  fFiducialEtaMin(-0.9),
  fFiducialEtaMax(0.9),
+ fFiducialPhiMin(0.),
+ fFiducialPhiMax(2*TMath::Pi()),
  fPtCut(2.0),
  fComment("No comment"),
  fDir(""),
index 2fbbf84..cbf8208 100755 (executable)
@@ -27,12 +27,15 @@ class AliJetReaderHeader : public TNamed
   virtual const char* GetBgDirectory(){return fBgDir.Data();}
   virtual const char* GetPattern() {return fPattern.Data();}
   virtual Float_t     GetFiducialEtaMin() const {return fFiducialEtaMin;}
-  virtual Float_t     GetFiducialEtaMax() const {return fFiducialEtaMax;}  
+  virtual Float_t     GetFiducialEtaMax() const {return fFiducialEtaMax;} 
+  virtual Float_t     GetFiducialPhiMin() const {return fFiducialPhiMin;}
+  virtual Float_t     GetFiducialPhiMax() const {return fFiducialPhiMax;}  
   virtual Float_t     GetPtCut() const {return fPtCut;}
   Int_t   GetNEvents()     const {return fLast-fFirst;}
   Int_t   GetFirstEvent()  const {return fFirst;}
   Int_t   GetLastEvent()   const {return fLast;}
   Int_t   GetDetector()    const {return fOption;}
+  Bool_t  GetDZ()          const {return fDZ;}
   Int_t   GetDebug()       const {return fDebug;}
   Int_t   GetSignalPerBg() const {return fSignalPerBg;}
          
@@ -46,18 +49,24 @@ class AliJetReaderHeader : public TNamed
   virtual void SetLastEvent(Int_t i=-1) {fLast=i;}
   virtual void SetFiducialEta(Float_t etamin, Float_t etamax) 
       { fFiducialEtaMin = etamin; fFiducialEtaMax = etamax;}
+  virtual void SetFiducialPhi(Float_t phimin, Float_t phimax) 
+      { fFiducialPhiMin = phimin; fFiducialPhiMax = phimax;}
   virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;}
-
+  virtual void SetDZ(Bool_t deadzone = 0) {fDZ = deadzone;}
   virtual void SetDetector(Int_t option = 0) {fOption = option;}
   virtual void SetDebug(Int_t debug = 0) {fDebug = debug;}
+
  protected:
   Int_t   fFirst;          // First and last events analyzed
   Int_t   fLast;           // in current set of files
   Int_t   fOption;         // detector used for jet reconstruction   
   Int_t   fDebug;          // debug option
+  Bool_t  fDZ;             // include dead zones or not 
   Int_t   fSignalPerBg;
   Float_t fFiducialEtaMin; // Fiducial minimum eta
   Float_t fFiducialEtaMax; // Fiducial maximum eta
+  Float_t fFiducialPhiMin; // Fiducial minimum phi
+  Float_t fFiducialPhiMax; // Fiducial maximum phi
   Float_t fPtCut;          // pt cut
   TString fComment;        // a comment
   TString fDir;            // directory with input files for signal
index 353a7e6..1c443d4 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-
-/* $Id$ */
-
-//_________________________________________________________________________
+//------------------------------------------------------------------
 //  Unit used by UA1 algorithm
-// --
-//*-- Author: Sarah Blyth (LBL/UCT)
-// --
-// Revised Version for JETAN - 30/03/2006
-// -- Magali Estienne (IReS)
+//  Authors: Sarah Blyth (LBL/UCT)
+//           Magali Estienne (IReS) (new version for JETAN)
+//------------------------------------------------------------------
 
 #include "AliJetUnitArray.h"
 
 
 ClassImp(AliJetUnitArray)
 
-
 AliJetUnitArray::AliJetUnitArray():
-    fUnitEnergy(0.),
-    fUnitEta(0.),
-    fUnitPhi(0.),
-    fUnitDeta(0.),
-    fUnitDphi(0.),
-    fUnitID(0),
-    fUnitNum(0),
-    fUnitClusterID(0),
-    fUnitFlag(kOutJet),
-    fUnitCutFlag(kPtSmaller),
-    fUnitSignalFlag(kGood),
-    fUnitDetectorFlag(kAll)
+  fUnitEnergy(0.0),
+  fUnitEta(0.0),
+  fUnitPhi(0.0),
+  fUnitDeta(0.),
+  fUnitDphi(0.),
+  fUnitID(0),
+  fUnitTrackID(0),
+  fUnitNum(0),
+  fUnitClusterID(0),
+  fUnitFlag(kOutJet),
+  fUnitCutFlag(kPtSmaller),
+  fUnitSignalFlag(kBad), 
+  fUnitDetectorFlag(kTpc),
+  fUnitPx(0.),
+  fUnitPy(0.),
+  fUnitPz(0.),
+  fUnitMass(0.)
 {
   // Default constructor
 }  
 
+AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t px, Float_t py, Float_t pz, Float_t Deta, Float_t Dphi, AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, AliJetFinderUnitCutFlagType_t cut, Float_t mass, Int_t clusId):
+  fUnitEnergy(en),
+  fUnitEta(eta),
+  fUnitPhi(phi),
+  fUnitDeta(Deta),
+  fUnitDphi(Dphi),
+  fUnitID(absId),
+  fUnitTrackID(esdId),
+  fUnitNum(0),
+  fUnitClusterID(clusId),
+  fUnitFlag(inout),
+  fUnitCutFlag(cut),
+  fUnitSignalFlag(kBad), 
+  fUnitDetectorFlag(det),
+  fUnitPx(px),
+  fUnitPy(py),
+  fUnitPz(pz),
+  fUnitMass(mass)
+{
+  // Constructor 2
+}
+
 AliJetUnitArray::~AliJetUnitArray()
 {
   // Destructor 
index 79d5fa5..20e810c 100644 (file)
@@ -1,19 +1,15 @@
-#ifndef ALIUA1JETFINDERUNIT_H
-#define ALIUA1JETFINDERUNIT_H
+#ifndef ALIJETUNITARRAY_H
+#define ALIJETUNITARRAY_H
 
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *  *  * See cxx source for full Copyright notice     */
 
-/* $Id$ */
-
-//======================================================================
-// Revised Version for JETAN - 08/11/2006
-// Functions added
-// Author: magali.estienne@ires.in2p3.fr
-//======================================================================
+//------------------------------------------------------------------
 //  Unit used by UA1 algorithm
-//
-//*-- Author: Sarah Blyth (LBL/UCT)
+//  Authors: Sarah Blyth (LBL/UCT)
+//           Magali Estienne (IReS) (new version for JETAN)
+//------------------------------------------------------------------
+
 
 #include <TObject.h>
 #include "AliJetFinderTypes.h"
@@ -22,6 +18,7 @@ class AliJetUnitArray : public TObject
 {
  public:
   AliJetUnitArray();
+  AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t px, Float_t py, Float_t pz, Float_t Deta, Float_t Dphi, AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, AliJetFinderUnitCutFlagType_t cut, Float_t mass, Int_t clusId);
   ~AliJetUnitArray();
 
   // Setter
@@ -31,6 +28,7 @@ class AliJetUnitArray : public TObject
   void SetUnitDeta(Float_t deta)      {fUnitDeta = deta;} 
   void SetUnitDphi(Float_t dphi)      {fUnitDphi = dphi;}
   void SetUnitID(Int_t id)            {fUnitID = id;}
+  void SetUnitTrackID(Int_t esdid)    {fUnitTrackID = esdid;}
   void SetUnitEntries(Int_t num)      {fUnitNum = num;}
   void SetUnitClusterID(Int_t id)     {fUnitClusterID = id;}
   void SetUnitFlag(AliJetFinderUnitFlagType_t flag)    
@@ -49,6 +47,8 @@ class AliJetUnitArray : public TObject
   { 
     fUnitDetectorFlag = detectorflag; 
   } 
+  void SetUnitPxPyPz(Double_t *pxyz)  {fUnitPx = pxyz[0]; fUnitPy = pxyz[1]; fUnitPz = pxyz[2];}
+  void SetUnitMass(Float_t mass) {fUnitMass = mass;}
 
   // Getter
   Float_t GetUnitEnergy() const            {return fUnitEnergy;}
@@ -57,8 +57,12 @@ class AliJetUnitArray : public TObject
   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(Double_t* p) const {p[0]=fUnitPx; p[1]=fUnitPy; p[2]=fUnitPz; return kTRUE;}
+
   AliJetFinderUnitFlagType_t GetUnitFlag() const     
   {
          return fUnitFlag;
@@ -76,25 +80,31 @@ class AliJetUnitArray : public TObject
          return fUnitDetectorFlag;
   }
 
+
   Bool_t operator>  ( AliJetUnitArray unit1) const;
   Bool_t operator<  ( AliJetUnitArray unit1) const;
   Bool_t operator== ( AliJetUnitArray unit1) const;
 
  protected:
-  Float_t         fUnitEnergy;        // Energy of the unit 
-  Float_t         fUnitEta;           // Eta of the unit
-  Float_t         fUnitPhi;           // Phi of the unit
-  Float_t         fUnitDeta;          // Delta Eta of the unit
-  Float_t         fUnitDphi;          // Delta Phi of the unit
-  Int_t           fUnitID;            // ID of the unit
-  Int_t           fUnitNum;           // number of units
-  Int_t           fUnitClusterID;     // ID of the unit
-  AliJetFinderUnitFlagType_t         fUnitFlag;          // Flag of the unit
-  AliJetFinderUnitCutFlagType_t      fUnitCutFlag;       // Flag of the unit
-  AliJetFinderUnitSignalFlagType_t   fUnitSignalFlag;    // Flag of the unit
-  AliJetFinderUnitDetectorFlagType_t fUnitDetectorFlag;  // Detector flag of the unit
+  Float_t         fUnitEnergy;                          // Energy of the unit 
+  Float_t         fUnitEta;                             // Eta of the unit
+  Float_t         fUnitPhi;                             // Phi of the unit
+  Float_t         fUnitDeta;                            // Delta Eta of the unit
+  Float_t         fUnitDphi;                            // Delta Phi of the unit
+  Int_t           fUnitID;                              // ID of the unit
+  Int_t           fUnitTrackID;                         // ID of the unit
+  Int_t           fUnitNum;                             // number of units
+  Int_t           fUnitClusterID;                       // ID of the unit
+  AliJetFinderUnitFlagType_t         fUnitFlag;         // Flag of the unit
+  AliJetFinderUnitCutFlagType_t      fUnitCutFlag;      // Flag of the unit
+  AliJetFinderUnitSignalFlagType_t   fUnitSignalFlag;   // Flag of the unit
+  AliJetFinderUnitDetectorFlagType_t fUnitDetectorFlag; // Detector flag of the unit
+  Float_t         fUnitPx;                              // Px of charged track
+  Float_t         fUnitPy;                              // Py of charged track
+  Float_t         fUnitPz;                              // Pz of charged track
+  Float_t         fUnitMass;                            // Mass of charged particle
 
-  ClassDef(AliJetUnitArray,3)
+  ClassDef(AliJetUnitArray,1)
 };
 
 #endif
diff --git a/JETAN/AliUA1JetFinderV2.cxx b/JETAN/AliUA1JetFinderV2.cxx
new file mode 100644 (file)
index 0000000..81da6b2
--- /dev/null
@@ -0,0 +1,878 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+/* $Id$ */
+//---------------------------------------------------------------------
+// UA1 Cone Algorithm Jet finder
+// manages the search for jets
+// Author: Rafael.Diaz.Valdes@cern.ch
+// (version in c++)
+// Modified to include neutral particles (magali.estienne@ires.in2p3.fr)
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+
+#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 "AliUA1JetFinderV2.h"
+#include "AliUA1JetHeaderV1.h"
+#include "AliJetUnitArray.h"
+#include "AliJetReaderHeader.h"
+#include "AliJetReader.h"
+#include "AliJet.h"
+#include "AliAODJet.h"
+
+
+ClassImp(AliUA1JetFinderV2)
+
+////////////////////////////////////////////////////////////////////////
+
+  AliUA1JetFinderV2::AliUA1JetFinderV2():
+    fDebug(0),
+    fOpt(0)
+{
+  // Constructor
+  fHeader = 0x0;
+  fLego   = 0x0;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+AliUA1JetFinderV2::~AliUA1JetFinderV2()
+
+{
+  // destructor
+}
+
+////////////////////////////////////////////////////////////////////////
+
+
+void AliUA1JetFinderV2::FindJets()
+
+{
+  //1) Fill cell map array
+  //2) calculate total energy and fluctuation level
+  //3) Run algorithm
+  //   3.1) look centroides in cell map
+  //   3.2) calculate total energy in cones
+  //   3.3) flag as a possible jet
+  //   3.4) reorder cones by energy
+  //4) subtract backg in accepted jets
+  //5) fill AliJet list
+
+  // transform input to pt,eta,phi plus lego
+    
+  AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+  TClonesArray*      fUnit = fReader->GetUnitArray();
+  Int_t              nCandidate = fReader->GetNumCandidate();
+  Int_t              nIn = fUnit->GetEntries();
+
+  if (nIn == 0) return;
+
+  // local arrays for input No Cuts
+  // Both pt < ptMin and pt > ptMin
+  Float_t* enT       = new Float_t[nCandidate];
+  Float_t* ptT       = new Float_t[nCandidate];
+  Float_t* etaT      = new Float_t[nCandidate];
+  Float_t* phiT      = new Float_t[nCandidate];
+  Float_t* detaT     = new Float_t[nCandidate];
+  Float_t* dphiT     = new Float_t[nCandidate];
+  Float_t* cFlagT    = new Float_t[nCandidate];
+  Float_t* cClusterT = new Float_t[nCandidate];
+  Float_t* idT       = new Float_t[nCandidate];
+  Int_t loop1 = 0;
+  Int_t* injet = new Int_t[nCandidate];
+  Int_t* sflag = new Int_t[nCandidate];
+
+
+  //total energy in array
+  Float_t  etbgTotal = 0.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
+  
+  // Information extracted from fUnitArray
+  // load input vectors and calculate total energy in array
+  for(Int_t i=0; i<nIn; i++) 
+    {
+      AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
+      if(uArray->GetUnitCutFlag()==1){ 
+       etCell[i] = uArray->GetUnitEnergy();
+       if (etCell[i] > 0.0) etCell[i] -= header->GetMinCellEt();
+       if (etCell[i] < 0.0) etCell[i] = 0.;
+       etaCell[i] = uArray->GetUnitEta();
+       phiCell[i] = uArray->GetUnitPhi();
+       flagCell[i] = 0; // default
+      }
+      else {
+       etCell[i]  = 0.;
+       etaCell[i] = uArray->GetUnitEta();
+       phiCell[i] = uArray->GetUnitPhi();
+       flagCell[i] = 0;
+      }
+
+      if(uArray->GetUnitEnergy()>0.){
+       ptT[loop1]   = uArray->GetUnitEnergy();
+       enT[loop1]   = uArray->GetUnitEnergy();
+       etaT[loop1]  = uArray->GetUnitEta();
+       phiT[loop1]  = uArray->GetUnitPhi();
+       detaT[loop1] = uArray->GetUnitDeta();
+       dphiT[loop1] = uArray->GetUnitDphi();
+       cFlagT[loop1]= uArray->GetUnitCutFlag();
+       idT[loop1]   = uArray->GetUnitID();
+       if(cFlagT[loop1] == 1) {
+         hPtTotal->Fill(ptT[loop1]);
+         //    fLego->Fill(etaT[i], phiT[i], ptT[i]);
+         etbgTotal+= ptT[loop1];
+       }       
+       loop1++;
+      }
+    }
+
+  //  fJets->SetNinput(nIn);
+  fJets->SetNinput(nCandidate);
+
+  // calculate total energy and fluctuation in map
+  Double_t meanpt = hPtTotal->GetMean();
+  Double_t ptRMS = hPtTotal->GetRMS();
+  Double_t npart = hPtTotal->GetEntries();
+  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* 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];
+  Int_t nJets; // to hold number of jets found by algorithm
+  Int_t nj;    // number of jets accepted
+  Float_t prec  = header->GetPrecBg();
+  Float_t bgprec = 1;
+  while(bgprec > prec){
+     //reset jet arrays in memory
+     memset(etaJet,0,sizeof(Float_t)*30);
+     memset(phiJet,0,sizeof(Float_t)*30);
+     memset(etJet,0,sizeof(Float_t)*30);
+     memset(etallJet,0,sizeof(Float_t)*30);
+     memset(etsigJet,0,sizeof(Float_t)*30);
+     memset(ncellsJet,0,sizeof(Int_t)*30);
+     memset(multJet,0,sizeof(Int_t)*30);
+     nJets = 0;
+     nj = 0;
+     // reset particles-jet array in memory
+     memset(injet,-1,sizeof(Int_t)*nCandidate);
+     //run cone algorithm finder
+     RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
+     //run background subtraction
+     if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
+       nj = header->GetNAcceptJets();
+     else
+       nj = nJets;
+     //subtract background
+     Float_t etbgTotalN = 0.0; //new background
+     if(header->GetBackgMode() == 1) // standar
+        SubtractBackg(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+     if(header->GetBackgMode() == 2) //cone
+        SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+     if(header->GetBackgMode() == 3) //ratio
+        SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+     if(header->GetBackgMode() == 4) //statistic
+        SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
+     //calc precision
+     if(etbgTotalN != 0.0)
+        bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
+     else
+        bgprec = 0;
+     etbgTotal = etbgTotalN; // update with new background estimation
+  } //end while
+
+  // add jets to list
+  Int_t* idxjets = new Int_t[nj];
+  Int_t nselectj = 0;
+  printf("Found %d jets \n", nj);
+  
+  for(Int_t kj=0; kj<nj; kj++){
+     if ((etaJet[kj] > (header->GetJetEtaMax())) ||
+          (etaJet[kj] < (header->GetJetEtaMin())) ||
+          (etJet[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
+      Float_t px, py,pz,en; // convert to 4-vector
+      px = etJet[kj] * TMath::Cos(phiJet[kj]);
+      py = etJet[kj] * TMath::Sin(phiJet[kj]);
+      pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj])));
+      en = TMath::Sqrt(px * px + py * py + pz * pz);
+      fJets->AddJet(px, py, pz, en);
+      AliAODJet jet(px, py, pz, en);
+      jet.Print("");
+      
+      AddJet(jet);
+      
+      idxjets[nselectj] = kj;
+      nselectj++;
+  }
+  //add signal percentage and total signal  in AliJets for analysis tool
+  Float_t* percentage  = new Float_t[nselectj];
+  Int_t* ncells      = new Int_t[nselectj];
+  Int_t* mult        = new Int_t[nselectj];
+  for(Int_t i = 0; i< nselectj; i++){
+     percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]];
+     ncells[i] = ncellsJet[idxjets[i]];
+     mult[i] = multJet[idxjets[i]];
+  }
+   //add particle-injet relationship ///
+  //   for(Int_t bj = 0; bj < nIn; bj++){
+   for(Int_t bj = 0; bj < nCandidate; bj++){
+       if(injet[bj] == -1) continue; //background particle
+       Int_t bflag = 0;
+       for(Int_t ci = 0; ci< nselectj; ci++){
+           if(injet[bj] == idxjets[ci]){
+              injet[bj]= ci;
+              bflag++;
+              break;
+           }
+       }
+       if(bflag == 0) injet[bj] = -1; // set as background particle
+   }
+  fJets->SetNCells(ncells);
+  fJets->SetPtFromSignal(percentage);
+  fJets->SetMultiplicities(mult);
+  fJets->SetInJet(injet);
+  fJets->SetEtaIn(etaT);
+  fJets->SetPhiIn(phiT);
+  fJets->SetPtIn(ptT);
+  fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi()));
+
+
+  //delete
+  delete enT;
+  delete ptT;
+  delete etaT;
+  delete phiT;
+  delete detaT;
+  delete dphiT;
+  delete cFlagT;
+  delete cClusterT;
+  delete idT;
+  delete injet;
+  delete sflag;
+  delete hPtTotal;
+  delete etCell;
+  delete etaCell;
+  delete phiCell;
+  delete flagCell;
+  delete etaJet;
+  delete phiJet;
+  delete etJet;
+  delete etsigJet;
+  delete etallJet;
+  delete ncellsJet;
+  delete multJet;
+  delete idxjets;
+  delete percentage;
+  delete ncells;
+  delete mult;
+
+
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, 
+                                   Int_t* flagCell, 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 nCell = nIn; 
+
+  AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+
+  // Parameters from header
+  Float_t minmove = header->GetMinMove();
+  Float_t maxmove = header->GetMaxMove();
+  Float_t rc      = header->GetRadius();
+  Float_t etseed  = header->GetEtSeed();
+
+  // tmp array of jets form algoritm
+  Float_t etaAlgoJet[30];
+  Float_t phiAlgoJet[30];
+  Float_t etAlgoJet[30];
+  Int_t   ncellsAlgoJet[30];
+
+  //run algorithm//
+
+  // sort cells by et
+  Int_t * index  = new Int_t[nCell];
+  TMath::Sort(nCell, etCell, index);
+
+  // variable used in centroide loop
+  Float_t eta = 0.0;
+  Float_t phi = 0.0;
+  Float_t eta0 = 0.0;
+  Float_t phi0 = 0.0;
+  Float_t etab = 0.0;
+  Float_t phib = 0.0;
+  Float_t etas = 0.0;
+  Float_t phis = 0.0;
+  Float_t ets = 0.0;
+  Float_t deta = 0.0;
+  Float_t dphi = 0.0;
+  Float_t dr = 0.0;
+  Float_t etsb = 0.0;
+  Float_t etasb = 0.0;
+  Float_t phisb = 0.0;
+
+
+  for(Int_t icell = 0; icell < nCell; icell++){
+        Int_t jcell = index[icell];
+        if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
+        if(flagCell[jcell] != 0) continue; // if cell was used before
+        eta  = etaCell[jcell];
+        phi  = phiCell[jcell];
+        eta0 = eta;
+        phi0 = phi;
+        etab = eta;
+        phib = phi;
+        ets  = etCell[jcell];
+        etas = 0.0;
+        phis = 0.0;
+        etsb = ets;
+        etasb = 0.0;
+        phisb = 0.0;
+        for(Int_t kcell =0; kcell < nCell; kcell++){
+            Int_t lcell = index[kcell];
+            if(lcell == jcell) continue; // cell itself
+            if(flagCell[lcell] != 0) continue; // cell used before
+            if(etCell[lcell] > etCell[jcell]) continue;
+            //calculate dr
+            deta = etaCell[lcell] - eta;
+                dphi = phiCell[lcell] - phi;
+                if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+                if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+                dr = TMath::Sqrt(deta * deta + dphi * dphi);
+            if(dr <= rc){
+               // calculate offset from initiate cell
+               deta = etaCell[lcell] - eta0;
+               dphi = phiCell[lcell] - phi0;
+               if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+                   if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+               etas = etas + etCell[lcell]*deta;
+               phis = phis + etCell[lcell]*dphi;
+               ets = ets + etCell[lcell];
+               //new weighted eta and phi including this cell
+               eta = eta0 + etas/ets;
+               phi = phi0 + phis/ets;
+               // if cone does not move much, just go to next step
+               dr = TMath::Sqrt((eta-etab)*(eta-etab) + (phi-phib)*(phi-phib));
+               if(dr <= minmove) break;
+               // cone should not move more than max_mov
+               dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
+               if(dr > maxmove){
+                  eta = etab;
+                  phi = phib;
+                  ets = etsb;
+                  etas = etasb;
+                  phis = phisb;
+               }else{ // store this loop information
+                 etab=eta;
+                 phib=phi;
+                 etsb = ets;
+                 etasb = etas;
+                 phisb = phis;
+               }
+            }
+        }//end of cells loop looking centroide
+
+        //avoid cones overloap (to be implemented in the future)
+
+        //flag cells in Rc, estimate total energy in cone
+        Float_t etCone = 0.0;
+        Int_t   nCellIn = 0;
+        rc = header->GetRadius();
+        for(Int_t ncell =0; ncell < nCell; ncell++){
+            if(flagCell[ncell] != 0) continue; // cell used before
+           //calculate dr
+            deta = etaCell[ncell] - eta;
+                dphi = phiCell[ncell] - phi;
+                if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+                if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+                dr = TMath::Sqrt(deta * deta + dphi * dphi);
+            if(dr <= rc){  // cell in cone
+               flagCell[ncell] = -1;
+               etCone+=etCell[ncell];
+               nCellIn++;
+            }
+        }
+
+        // select jets with et > background
+        // estimate max fluctuation of background in cone
+        Double_t ncellin = (Double_t)nCellIn;
+        Double_t ntcell  = (Double_t)nCell;
+        Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
+        // min cone et
+        Double_t etcmin = etCone ;  // could be used etCone - etmin !!
+        //desicions !! etbmax < etcmin
+        for(Int_t mcell =0; mcell < nCell; mcell++){
+            if(flagCell[mcell] == -1){
+              if(etbmax < etcmin)
+                 flagCell[mcell] = 1; //flag cell as used
+              else
+                 flagCell[mcell] = 0; // leave it free
+            }
+        }
+        //store tmp jet info !!!
+       if(etbmax < etcmin) {
+             etaAlgoJet[nJets] = eta;
+             phiAlgoJet[nJets] = phi;
+             etAlgoJet[nJets] = etCone;
+             ncellsAlgoJet[nJets] = nCellIn;
+             nJets++;
+        }
+
+  } // end of cells loop
+
+  //reorder jets by et in cone
+  //sort jets by energy
+  Int_t * idx  = new Int_t[nJets];
+  TMath::Sort(nJets, etAlgoJet, idx);
+  for(Int_t p = 0; p < nJets; p++){
+     etaJet[p] = etaAlgoJet[idx[p]];
+     phiJet[p] = phiAlgoJet[idx[p]];
+     etJet[p] = etAlgoJet[idx[p]];
+     etallJet[p] = etAlgoJet[idx[p]];
+     ncellsJet[p] = ncellsAlgoJet[idx[p]];
+  }
+
+
+  //delete
+  delete index;
+  delete idx;
+
+}
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinderV2::SubtractBackg(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,
+                      Int_t* multJet, Int_t* injet)
+{
+  //background subtraction using cone method but without correction in dE/deta distribution
+
+  //calculate energy inside and outside cones
+  AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+  Float_t rc= header->GetRadius();
+  Float_t etIn[30];
+  Float_t etOut = 0;
+  for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+     // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
+     for(Int_t ijet=0; ijet<nJ; ijet++){
+         Float_t deta = etaT[jpart] - etaJet[ijet];
+             Float_t dphi = phiT[jpart] - phiJet[ijet];
+         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+             if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+             Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+         if(dr <= rc){ // particles inside this cone
+             multJet[ijet]++;
+             injet[jpart] = ijet;
+             if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
+                etIn[ijet] += ptT[jpart];
+                if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
+             }
+             break;
+         }
+     }// end jets loop
+     if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
+        etOut += ptT[jpart]; // particle outside cones and pt cut
+  } //end particle loop
+
+  //estimate jets and background areas
+  Float_t areaJet[30];
+  Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
+  for(Int_t k=0; k<nJ; k++){
+      Float_t detamax = etaJet[k] + rc;
+      Float_t detamin = etaJet[k] - rc;
+      Float_t accmax = 0.0; Float_t accmin = 0.0;
+      if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
+         Float_t h = header->GetLegoEtaMax() - etaJet[k];
+         accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+      }
+      if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
+         Float_t h = header->GetLegoEtaMax() + etaJet[k];
+         accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+      }
+      areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
+      areaOut = areaOut - areaJet[k];
+  }
+  //subtract background using area method
+  for(Int_t ljet=0; ljet<nJ; ljet++){
+     Float_t areaRatio = areaJet[ljet]/areaOut;
+     etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
+  }
+
+  // estimate new total background
+  Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
+  etbgTotalN = etOut*areaT/areaOut;
+
+
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinderV2::SubtractBackgStat(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,
+                      Int_t* multJet, Int_t* injet)
+{
+
+  //background subtraction using statistical method
+  AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+  Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
+
+  //calculate energy inside
+  Float_t rc= header->GetRadius();
+  Float_t etIn[30];
+
+  for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+     //if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
+     for(Int_t ijet=0; ijet<nJ; ijet++){
+         Float_t deta = etaT[jpart] - etaJet[ijet];
+             Float_t dphi = phiT[jpart] - phiJet[ijet];
+         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+             if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+             Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+         if(dr <= rc){ // particles inside this cone
+             multJet[ijet]++;
+             injet[jpart] = ijet;
+             if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
+                etIn[ijet]+= ptT[jpart];
+                if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
+             }
+             break;
+         }
+     }// end jets loop
+  } //end particle loop
+
+  //calc jets areas
+  Float_t areaJet[30];
+  Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
+  for(Int_t k=0; k<nJ; k++){
+      Float_t detamax = etaJet[k] + rc;
+      Float_t detamin = etaJet[k] - rc;
+      Float_t accmax = 0.0; Float_t accmin = 0.0;
+      if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
+         Float_t h = header->GetLegoEtaMax() - etaJet[k];
+         accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+      }
+      if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
+         Float_t h = header->GetLegoEtaMax() + etaJet[k];
+         accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
+      }
+      areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
+  }
+
+  //subtract background using area method
+  for(Int_t ljet=0; ljet<nJ; ljet++){
+     Float_t areaRatio = areaJet[ljet]/areaOut;
+     etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
+  }
+
+  etbgTotalN = etbgStat;
+
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinderV2::SubtractBackgCone(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,
+                      Int_t* multJet, Int_t* injet)
+{
+   // Cone background subtraction method taking into acount dEt/deta distribution
+    AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+   //general
+   Float_t rc= header->GetRadius();
+   Float_t etamax = header->GetLegoEtaMax();
+   Float_t etamin = header->GetLegoEtaMin();
+   Int_t ndiv = 100;
+
+   // jet energy and area arrays
+   TH1F* hEtJet[30];
+   TH1F* hAreaJet[30];
+   for(Int_t mjet=0; mjet<nJ; mjet++){
+     char hEtname[256]; char hAreaname[256];
+     sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
+     hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
+     hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
+  }
+   // background energy and area
+   TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
+   TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
+
+   //fill energies
+   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
+         Float_t deta = etaT[jpart] - etaJet[ijet];
+        Float_t dphi = phiT[jpart] - phiJet[ijet];
+         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+        if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+        Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+         if(dr <= rc){ // particles inside this cone
+             injet[jpart] = ijet;
+             multJet[ijet]++;
+             if((fReader->GetCutFlag(jpart)) == 1){// pt cut
+                hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
+                if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
+             }
+             break;
+         }
+     }// end jets loop
+     if(injet[jpart] == -1  && fReader->GetCutFlag(jpart) == 1)
+        hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+  } //end particle loop
+
+   //calc areas
+   Float_t eta0 = etamin;
+   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
+   Float_t eta1 = eta0 + etaw;
+   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
+      Float_t etac = eta0 + etaw/2.0;
+      Float_t areabg = etaw*2.0*TMath::Pi();
+      for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
+          Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
+          Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
+          Float_t acc0 = 0.0; Float_t acc1 = 0.0;
+          Float_t areaj = 0.0;
+          if(deta0 > rc && deta1 < rc){
+             acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+             areaj = acc1;
+          }
+          if(deta0 < rc && deta1 > rc){
+             acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+             areaj = acc0;
+          }
+          if(deta0 < rc && deta1 < rc){
+             acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+             acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+             if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
+             if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
+             if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
+          }
+          hAreaJet[ijet]->Fill(etac,areaj);
+          areabg = areabg - areaj;
+      } // end jets loop
+      hAreaBackg->Fill(etac,areabg);
+      eta0 = eta1;
+      eta1 = eta1 + etaw;
+   } // end loop for all eta bins
+
+   //subtract background
+   for(Int_t kjet=0; kjet<nJ; kjet++){
+       etJet[kjet] = 0.0; // first  clear etJet for this jet
+       for(Int_t bin = 0; bin< ndiv; bin++){
+           if(hAreaJet[kjet]->GetBinContent(bin)){
+              Float_t areab = hAreaBackg->GetBinContent(bin);
+              Float_t etb = hEtBackg->GetBinContent(bin);
+              Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
+              etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
+           }
+       }
+   }
+
+   // calc background total
+   Double_t etOut = hEtBackg->Integral();
+   Double_t areaOut = hAreaBackg->Integral();
+   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
+   etbgTotalN = etOut*areaT/areaOut;
+
+   //delete
+   for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
+       delete hEtJet[ljet];
+       delete hAreaJet[ljet];
+   }
+
+   delete hEtBackg;
+   delete hAreaBackg;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+
+void AliUA1JetFinderV2::SubtractBackgRatio(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,
+                       Int_t* multJet, Int_t* injet)
+{
+   // Ratio background subtraction method taking into acount dEt/deta distribution
+    AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+   //factor F calc before
+    Float_t bgRatioCut = header->GetBackgCutRatio();
+
+
+   //general
+   Float_t rc= header->GetRadius();
+   Float_t etamax = header->GetLegoEtaMax();
+   Float_t etamin = header->GetLegoEtaMin();
+   Int_t ndiv = 100;
+
+   // jet energy and area arrays
+   TH1F* hEtJet[30];
+   TH1F* hAreaJet[30];
+   for(Int_t mjet=0; mjet<nJ; mjet++){
+     char hEtname[256]; char hAreaname[256];
+     sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
+     hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);        // change range
+     hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);  // change range
+  }
+   // background energy and area
+   TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);         // change range
+   TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);  // change range
+
+   //fill energies
+   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
+     //if((fReader->GetCutFlag(jpart)) != 1) continue;
+     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
+         Float_t deta = etaT[jpart] - etaJet[ijet];
+             Float_t dphi = phiT[jpart] - phiJet[ijet];
+         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+             if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+             Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+         if(dr <= rc){ // particles inside this cone
+            multJet[ijet]++;
+            injet[jpart] = ijet;
+            if((fReader->GetCutFlag(jpart)) == 1){ //pt cut
+               hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
+               if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
+            }
+            break;
+         }
+     }// end jets loop
+     if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
+  } //end particle loop
+
+   //calc areas
+   Float_t eta0 = etamin;
+   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
+   Float_t eta1 = eta0 + etaw;
+   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
+      Float_t etac = eta0 + etaw/2.0;
+      Float_t areabg = etaw*2.0*TMath::Pi();
+      for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
+          Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
+          Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
+          Float_t acc0 = 0.0; Float_t acc1 = 0.0;
+          Float_t areaj = 0.0;
+          if(deta0 > rc && deta1 < rc){
+             acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+             areaj = acc1;
+          }
+          if(deta0 < rc && deta1 > rc){
+             acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+             areaj = acc0;
+          }
+          if(deta0 < rc && deta1 < rc){
+             acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
+             acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
+             if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
+             if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
+             if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
+          }
+          hAreaJet[ijet]->Fill(etac,areaj);
+          areabg = areabg - areaj;
+      } // end jets loop
+      hAreaBackg->Fill(etac,areabg);
+      eta0 = eta1;
+      eta1 = eta1 + etaw;
+   } // end loop for all eta bins
+
+   //subtract background
+   for(Int_t kjet=0; kjet<nJ; kjet++){
+       etJet[kjet] = 0.0; // first  clear etJet for this jet
+       for(Int_t bin = 0; bin< ndiv; bin++){
+           if(hAreaJet[kjet]->GetBinContent(bin)){
+              Float_t areab = hAreaBackg->GetBinContent(bin);
+              Float_t etb = hEtBackg->GetBinContent(bin);
+              Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
+              etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
+           }
+       }
+   }
+
+   // calc background total
+   Double_t etOut = hEtBackg->Integral();
+   Double_t areaOut = hAreaBackg->Integral();
+   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
+   etbgTotalN = etOut*areaT/areaOut;
+
+   //delete
+   for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
+       delete hEtJet[ljet];
+       delete hAreaJet[ljet];
+   }
+
+   delete hEtBackg;
+   delete hAreaBackg;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+
+void AliUA1JetFinderV2::Reset()
+{
+  fLego->Reset();
+  fJets->ClearJets();
+  AliJetFinder::Reset();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinderV2::WriteJHeaderToFile()
+{
+  AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+  header->Write();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliUA1JetFinderV2::Init()
+{
+  // initializes some variables
+  AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
+   // book lego
+  fLego = new
+    TH2F("legoH","eta-phi",
+        header->GetLegoNbinEta(), header->GetLegoEtaMin(),
+        header->GetLegoEtaMax(),  header->GetLegoNbinPhi(),
+        header->GetLegoPhiMin(),  header->GetLegoPhiMax());
+
+  fDebug = fReader->GetReaderHeader()->GetDebug();
+  fOpt = fReader->GetReaderHeader()->GetDetector();
+
+  if(fOpt>0)
+    fReader->CreateTasks();
+
+}
diff --git a/JETAN/AliUA1JetFinderV2.h b/JETAN/AliUA1JetFinderV2.h
new file mode 100644 (file)
index 0000000..b382345
--- /dev/null
@@ -0,0 +1,69 @@
+#ifndef ALIUA1JETFINDERV2_H
+#define ALIUA1JETFINDERV2_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+//---------------------------------------------------------------------
+// UA1 Cone Algorithm Finder V1
+// manages the search for jets
+// Author: Rafael.Diaz.Valdes@cern.ch
+// (version in c++)
+// Modified to include neutral particles (magali.estienne@ires.in2p3.fr)
+//---------------------------------------------------------------------
+
+#include "AliJetFinder.h"
+class AliUA1JetHeaderV1;
+class TH2F;
+
+class AliUA1JetFinderV2 : public AliJetFinder
+{
+ public:
+
+  AliUA1JetFinderV2();
+  ~AliUA1JetFinderV2();
+
+  // others
+  void FindJets();
+
+  void RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, 
+                  Int_t* flagCell, 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 SubtractBackg(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,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* etJet,Float_t* etaJet, 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* etJet,Float_t* etaJet, 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* etJet,Float_t* etaJet, Float_t* phiJet,
+                      Float_t* etsigJet, Int_t* multJet, Int_t* injet);
+  void Reset();
+  void Init();
+  void WriteJHeaderToFile();
+
+ protected:
+  AliUA1JetFinderV2(const AliUA1JetFinderV2& rJetF1);
+  AliUA1JetFinderV2& operator = (const AliUA1JetFinderV2& rhsf);
+  TH2F           * fLego;           //Lego Histo
+  Int_t fDebug;
+  Int_t fOpt;
+
+  ClassDef(AliUA1JetFinderV2,1)
+};
+
+#endif
index 3933bfc..73b1f15 100644 (file)
@@ -35,6 +35,7 @@ AliUA1JetHeaderV1::AliUA1JetHeaderV1():
     fConeRadius(0.3),
     fEtSeed(3.0),
     fMinJetEt(10.),
+    fMinCellEt(0.0), // Temporarily added -> To be removed if not necessary
     fMinMove(0.05),
     fMaxMove(0.15),
     fBackgMode(1),
index 69f2adf..f54eba9 100644 (file)
@@ -26,6 +26,7 @@ class AliUA1JetHeaderV1 : public AliJetHeader
   Float_t GetMaxMove()   const  {return fMaxMove;}
   Float_t GetEtSeed()    const  {return fEtSeed;}
   Float_t GetMinJetEt()  const  {return fMinJetEt;}
+  Float_t GetMinCellEt() const  {return fMinCellEt;} // Added temporarily !!! To be removed if not necessary
   Int_t   GetLegoNbinEta()   const  {return fLegoNbinEta;}
   Int_t   GetLegoNbinPhi()   const  {return fLegoNbinPhi;}
   Float_t GetLegoEtaMin()    const  {return fLegoEtaMin;}
@@ -45,6 +46,7 @@ class AliUA1JetHeaderV1 : public AliJetHeader
   void SetMaxMove(Float_t f) {fMaxMove=f;}
   void SetEtSeed(Float_t f) {fEtSeed=f;}
   void SetMinJetEt(Float_t f) {fMinJetEt=f;}
+  void SetMinCellEt(Float_t f) {fMinCellEt=f;}
   void SetLegoNbinEta(Int_t f) {fLegoNbinEta=f;}
   void SetLegoNbinPhi(Int_t f) {fLegoNbinPhi=f;}
   void SetLegoEtaMin(Float_t f) {fLegoEtaMin=f;}
@@ -66,6 +68,7 @@ protected:
   Float_t fConeRadius;      //  Cone radius
   Float_t fEtSeed;          //  Min. Et for seed
   Float_t fMinJetEt;        //  Min Et of jet
+  Float_t fMinCellEt;       //  Min Et in one cell
   // parameters of backgound substraction
   Float_t fMinMove;         // min cone move
   Float_t fMaxMove;         // max cone move
index e4668bb..0731d4f 100644 (file)
@@ -19,6 +19,7 @@
 #pragma link C++ class AliJetAnalysis+;
 #pragma link C++ class AliJetDistributions+;
 #pragma link C++ class AliUA1JetFinderV1+;
+#pragma link C++ class AliUA1JetFinderV2+;
 #pragma link C++ class AliUA1JetHeaderV1+;
 #pragma link C++ class AliJetGrid+;
 #pragma link C++ class AliJetUnitArray+;
@@ -26,6 +27,7 @@
 #pragma link C++ class AliJetHadronCorrectionv0+;
 #pragma link C++ class AliJetHadronCorrectionv1+; 
 #pragma link C++ class AliJetFillUnitArrayTracks+;
+#pragma link C++ class AliJetFillUnitArrayEMCalDigits+;
 #pragma link C++ class AliJetDummyGeo+;
 #pragma link C++ class AliAnalysisTaskJets+;
 #pragma link C++ class AliDAJetFinder+;
index 1fe5d25..9810012 100644 (file)
@@ -7,10 +7,10 @@ SRCS =        AliJet.cxx AliJetHeader.cxx \
        AliLeading.cxx \
        AliJetProductionData.cxx AliJetProductionDataPDC2004.cxx \
        AliJetAnalysis.cxx AliJetDistributions.cxx \
-       AliUA1JetFinderV1.cxx AliUA1JetHeaderV1.cxx \
+       AliUA1JetFinderV1.cxx AliUA1JetFinderV2.cxx AliUA1JetHeaderV1.cxx \
        AliJetGrid.cxx AliJetUnitArray.cxx AliJetHadronCorrection.cxx \
        AliJetHadronCorrectionv0.cxx AliJetHadronCorrectionv1.cxx \
-       AliJetFillUnitArrayTracks.cxx \
+       AliJetFillUnitArrayTracks.cxx AliJetFillUnitArrayEMCalDigits.cxx\
        AliJetDummyGeo.cxx \
        AliJetFinderTypes.cxx \
        AliAnalysisTaskJets.cxx \