- Clean-up of data members and methods in base classes and derived classes
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 19 Dec 2006 13:00:55 +0000 (13:00 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 19 Dec 2006 13:00:55 +0000 (13:00 +0000)
- Use of UnitArray to read TPC and EMCAL information. (Magali Estienne)

24 files changed:
JETAN/AliJetESDReader.cxx
JETAN/AliJetESDReader.h
JETAN/AliJetESDReaderHeader.cxx
JETAN/AliJetESDReaderHeader.h
JETAN/AliJetESDmcReader.cxx
JETAN/AliJetESDmcReader.h
JETAN/AliJetFillUnitArrayTracks.cxx
JETAN/AliJetFillUnitArrayTracks.h
JETAN/AliJetFinder.cxx
JETAN/AliJetFinder.h
JETAN/AliJetKineReader.cxx
JETAN/AliJetKineReader.h
JETAN/AliJetMCReader.cxx
JETAN/AliJetMCReader.h
JETAN/AliJetMCReaderHeader.cxx
JETAN/AliJetMCReaderHeader.h
JETAN/AliJetReader.cxx
JETAN/AliJetReader.h
JETAN/AliJetReaderHeader.cxx
JETAN/AliJetReaderHeader.h
JETAN/AliUA1JetFinder.cxx
JETAN/AliUA1JetFinder.h
JETAN/JETANLinkDef.h
JETAN/libJETAN.pkg

index 9c2a988..364c8b3 100755 (executable)
@@ -16,7 +16,8 @@
 //------------------------------------------------------------------------- 
 // Jet ESD Reader 
 // ESD reader for jet analysis
-// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+// Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+//          Magali Estienne <magali.estienne@IReS.in2p3.fr>
 //------------------------------------------------------------------------- 
 
 
 #include <TSystem.h>
 #include <TLorentzVector.h>
 #include <TVector3.h>
+#include <TGeoManager.h>
+
 #include "AliJetESDReader.h"
 #include "AliJetESDReaderHeader.h"
 #include "AliESD.h"
 #include "AliESDtrack.h"
+#include "AliEMCALGeometry.h"
+#include "AliJetFillUnitArrayTracks.h"
+#include "AliJetUnitArray.h"
 
 ClassImp(AliJetESDReader)
 
 AliJetESDReader::AliJetESDReader():
-  fMass(0),
-  fSign(0)
+    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) 
 {
   // Constructor    
-  fReaderHeader = 0x0;
+    SetEMCALGeometry();
 }
 
 //____________________________________________________________________________
@@ -44,6 +64,10 @@ AliJetESDReader::AliJetESDReader():
 AliJetESDReader::~AliJetESDReader()
 {
   // Destructor
+    delete fChain;
+    delete fESD;
+    delete fTpcGrid;
+    delete fEmcalGrid;
 }
 
 //____________________________________________________________________________
@@ -61,7 +85,7 @@ void AliJetESDReader::OpenInputFiles()
   
    void *dir  = gSystem->OpenDirectory(dirName);
    const char *name = 0x0;
-   int nesd = fReaderHeader->GetNesd();
+   int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
    int a = 0;
    while ((name = gSystem->GetDirEntry(dir))){
        if (a>=nesd) continue;
@@ -94,6 +118,7 @@ void AliJetESDReader::OpenInputFiles()
 }
 
 void AliJetESDReader::ConnectTree(TTree* tree) {
+    // Connect the tree
      fChain = (TChain*) tree;
      
      fChain->SetBranchAddress("ESD",    &fESD);
@@ -121,7 +146,8 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t event)
   
   // clear momentum array
    ClearArray();
-  
+   fDebug = fReaderHeader->GetDebug();
+   InitParameters();
   // get event from chain
   fChain->GetTree()->GetEntry(event);
 
@@ -149,14 +175,14 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t event)
       track->GetPxPyPz(mom);
       p3.SetXYZ(mom[0],mom[1],mom[2]);
       pt = p3.Pt();
-      if ((status & AliESDtrack::kTPCrefit) == 0) continue;    // quality check
-      if ((status & AliESDtrack::kITSrefit) == 0) continue;    // quality check
+      if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
+      if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() 
-         && TMath::Abs(track->GetLabel()) > 10000)  continue;   // quality check
+         && TMath::Abs(track->GetLabel()) > 10000)  continue;      // quality check
       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() 
-         && TMath::Abs(track->GetLabel()) < 10000)  continue;   // quality check
+         && TMath::Abs(track->GetLabel()) < 10000)  continue;      // quality check
       eta = p3.Eta();
-      if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
+      if ( (eta > etaMax) || (eta < etaMin))         continue;      // checking eta cut
       
       new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
       sflag[goodTrack]=0;
@@ -168,9 +194,64 @@ Bool_t AliJetESDReader::FillMomentumArray(Int_t event)
   // set the signal flags
   fSignalFlag.Set(goodTrack,sflag);
   fCutFlag.Set(goodTrack,cflag);
+
+//
+//
+  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");
+  }
+  delete fillUAFromTracks;
   return kTRUE;
 }
 
+
+void AliJetESDReader::SetEMCALGeometry()
+{
+  // Define EMCAL geometry to be able to read ESDs
+    fGeom = AliEMCALGeometry::GetInstance();
+    if (fGeom == 0)
+       fGeom = AliEMCALGeometry::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
+    
+    // To be setted to run some AliEMCALGeometry functions
+    TGeoManager::Import("geometry.root");
+    fGeom->GetTransformationForSM();  
+    
+    printf("\n EMCal Geometry set ! \n");
+
+}
   
+void AliJetESDReader::InitParameters()
+{
+    // Initialise parameters
+    fHCorrection    = 0;                 // For hadron correction
+    fHadCorr        = 0;                 // For hadron correction
+    fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
+    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;
+}
 
 
index 790d991..e3d8f58 100755 (executable)
@@ -8,30 +8,59 @@
 // Jet ESD Reader 
 // ESD reader for jet analysis
 // Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+//=========================================================================
+// Modified in order to use a fUnitArray object instead of a fMomentumArray
+// Includes EMCal Geometry, fUnitArray, grid objects and tools for Hadron correction
+// Author : magali.estienne@ires.in2p3.fr
 //---------------------------------------------------------------------
 
 #include "AliJetReader.h"
+#include "AliJetUnitArray.h"
+#include "AliJetGrid.h"
 class AliJetESDReaderHeader;
-
+class AliEMCALGeometry;
+class AliJetHadronCorrection;
+class AliJetUnitArray;
+class AliJetReaderHeader;
+class AliESD;
 
 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   ConnectTree(TTree* tree);
-  
+  virtual void SetTPCGrid(AliJetGrid *grid)   {fTpcGrid = grid;}
+  virtual void SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;}
+  // Correction of hadronic energy
+  virtual void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;}
+  virtual void SetHadronCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;}
+ private:
+  void SetEMCALGeometry();
+  void InitParameters();
  protected:
-  Float_t fMass;    // Particle mass
-  Int_t   fSign;    // Particle sign
+  AliEMCALGeometry           *fGeom;             //!EMCAL Geometry 
+  TChain                     *fChain;            // chain for reconstructed tracks
+  AliESD                     *fESD;              // pointer to esd
+  AliJetHadronCorrectionv1   *fHadCorr;          // Pointer to Hadron Correction Object 
+  AliJetGrid                 *fTpcGrid;          // Pointer to grid object
+  AliJetGrid                 *fEmcalGrid;        // 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
+  Int_t                       fNIn;              // Number of Array filled in UnitArray
+  Int_t                       fOpt;              // Detector to be used for jet reconstruction
+  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)
 };
index 22ddb93..67e0c23 100755 (executable)
@@ -27,7 +27,8 @@ ClassImp(AliJetESDReaderHeader)
 AliJetESDReaderHeader::AliJetESDReaderHeader():
   AliJetReaderHeader("AliJetESDReaderHeader"), 
   fReadSignalOnly(kFALSE),
-  fReadBkgdOnly(kFALSE)
+  fReadBkgdOnly(kFALSE),
+  fNesd(0)
 {
   // Default constructor
   
index a0c36b8..426b504 100755 (executable)
@@ -18,17 +18,27 @@ class AliJetESDReaderHeader : public AliJetReaderHeader
   virtual ~AliJetESDReaderHeader();
   
   // Getters
+  Float_t GetDCA()        const {return fDCA;}       
+  Float_t GetTLength()    const {return fTLength;}    
   Bool_t  ReadSignalOnly() const  {return fReadSignalOnly;}
   Bool_t  ReadBkgdOnly() const  {return fReadBkgdOnly;}
-         
+  Int_t   GetNesd()       const {return fNesd;}
+           
   // Setters
+  virtual void SetDCA(Float_t dca = 0.0) {fDCA = dca;}
+  virtual void SetTLength(Float_t length = 0.0) {fTLength = length;}
   virtual void SetReadSignalOnly(Bool_t flag = kTRUE) {fReadSignalOnly = flag;}
   virtual void SetReadBkgdOnly(Bool_t flag = kTRUE) {fReadBkgdOnly = flag;}
-  
+  virtual void SetNumberOfESD(Int_t i=1) {fNesd = i;}
+    
  protected:
   //parameters set by user
+  Float_t fDCA;            // dca cut
+  Float_t fTLength;        // track length cut
   Bool_t  fReadSignalOnly; // read particles from signal event only
   Bool_t  fReadBkgdOnly;   // read particles from bkgd event only
+  Int_t   fNesd;           // number of esds
+  
   ClassDef(AliJetESDReaderHeader,2);
 };
  
index ad10e5d..5b6318e 100644 (file)
 ClassImp(AliJetESDmcReader)
 
 AliJetESDmcReader::AliJetESDmcReader():
-    fAliHeader(0),
-    fMass(0),
-    fSign(0)
+    AliJetESDReader(),
+    fChainMC(0x0),
+    fAliHeader(0)
 {
   // Constructor    
-  fReaderHeader = 0x0;
 }
 
 //____________________________________________________________________________
@@ -47,6 +46,7 @@ AliJetESDmcReader::AliJetESDmcReader():
 AliJetESDmcReader::~AliJetESDmcReader()
 {
   // Destructor
+    delete fChainMC;
 }
 
 //____________________________________________________________________________
@@ -65,7 +65,7 @@ void AliJetESDmcReader::OpenInputFiles()
   // Add files matching patters to the chain
   void *dir  = gSystem->OpenDirectory(dirName);
   const char *name = 0x0;
-  int nesd = fReaderHeader->GetNesd();
+  int nesd = ((AliJetESDReaderHeader*)fReaderHeader)->GetNesd();
   int a = 0;
   while ((name = gSystem->GetDirEntry(dir))){
     if (a>=nesd) continue;
index 83662df..03f9ef2 100644 (file)
@@ -8,31 +8,24 @@
 // ESD reader for jet analysis (it reads the esd and the MC trees)
 // Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
 
-#include "AliJetReader.h"
+#include "AliJetESDReader.h"
 class AliJetESDReaderHeader;
 class AliHeader;
 
 
 
-class AliJetESDmcReader : public AliJetReader
+class AliJetESDmcReader : public AliJetESDReader
 {
  public: 
   AliJetESDmcReader();
   virtual ~AliJetESDmcReader();
-
-  // 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();
    
  protected:
+  TChain*    fChainMC;   // chain for mc information
   AliHeader* fAliHeader; //! Event header
-  Float_t    fMass;      //! Particle mass
-  Int_t      fSign;      //! Particle sign
-
   ClassDef(AliJetESDmcReader,1)
 };
  
index b41ddac..d1f059c 100644 (file)
@@ -80,7 +80,7 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
   : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information")
 {
   // constructor
-   fNIn = 0;
+  fNIn = 0;
   fOpt = 0;
   fDebug = 0;
   fNphi = 0;
@@ -99,11 +99,11 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
 //____________________________________________________________________________
 void AliJetFillUnitArrayTracks::SetEMCALGeometry()
 {
-  // Set EMCAL geometry information
-  fGeom = AliEMCALGeometry::GetInstance();
-  if (fGeom == 0)
-    fGeom = AliEMCALGeometry::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
-  if(fDebug>1) printf("\n EMCAL Geometry setted ! \n");
+    // Set EMCAL geometry information
+    fGeom = AliEMCALGeometry::GetInstance();
+    if (fGeom == 0)
+       fGeom = AliEMCALGeometry::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
+    if(fDebug>1) printf("\n EMCAL Geometry setted ! \n");
 }
 
 //____________________________________________________________________________
index 4dd0c6a..f0d1697 100644 (file)
@@ -59,7 +59,7 @@ 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                       fHCorrection;     // Hadron correction flag
   Int_t                       fNIn;             // Number of Array filled in UnitArray
   Int_t                       fOpt;             // Detector to be used for jet reconstruction
   Int_t                       fDebug;           // Debug option
index 85316aa..156aabc 100644 (file)
@@ -135,21 +135,35 @@ void AliJetFinder::Run()
   // write headers
   WriteHeaders();
   // loop over events
-  Int_t nFirst,nLast;
+  Int_t nFirst, nLast, option, debug, arrayInitialised; 
   nFirst = fReader->GetReaderHeader()->GetFirstEvent();
-  nLast = fReader->GetReaderHeader()->GetLastEvent();
+  nLast  = fReader->GetReaderHeader()->GetLastEvent();
+
+  option = fReader->GetReaderHeader()->GetDetector();
+  debug  = fReader->GetReaderHeader()->GetDebug();
+  arrayInitialised = fReader->GetArrayInitialised();
+
   // loop over events
   for (Int_t i=nFirst;i<nLast;i++) {
       fReader->FillMomentumArray(i);
       fLeading->FindLeading(fReader);
       fReader->GetGenJets(fGenJets);
-      FindJets();
+
+      if (option == 0) { // TPC with fMomentumArray
+         if(debug > 1) 
+             printf("In FindJetsTPC() routine: find jets with fMomentumArray !!!\n");
+         FindJetsTPC();
+      } else {
+          if(debug > 1) printf("In FindJets() routine: find jets with fUnitArray !!!\n");
+          FindJets();
+      }
       if (fOut) WriteJetsToFile(i);
       if (fPlots) fPlots->FillHistos(fJets);
       fLeading->Reset();
       fGenJets->ClearJets();
       Reset();
   } 
+
   // write out
   if (fPlots) {
       fPlots->Normalize();
index dbb18f5..a73091a 100755 (executable)
@@ -44,8 +44,9 @@ class AliJetFinder : public TObject
   virtual void   WriteRHeaderToFile();  
   // the following have to be implemented for each specific finder
   virtual void Init() {}
-  virtual void Reset() { }
-  virtual void FindJets() { }
+  virtual void Reset() {}
+  virtual void FindJets() {}
+  virtual void FindJetsTPC(){}
   virtual void WriteJHeaderToFile() { }
   // some methods to allow steering from the outside
   virtual Bool_t ProcessEvent(Long64_t entry);
index 7797e8e..55bf580 100644 (file)
 ClassImp(AliJetKineReader)
 
 AliJetKineReader::AliJetKineReader():
-  fRunLoader(0x0),
-  fAliHeader(0x0),
-  fMass(0),
-  fPdgC(0)
+    AliJetReader(),  
+    fRunLoader(0x0),
+    fAliHeader(0x0)
 {
   // Default constructor
-  fReaderHeader = 0x0;
 }
 
 //____________________________________________________________________________
@@ -55,7 +53,6 @@ AliJetKineReader::AliJetKineReader():
 AliJetKineReader::~AliJetKineReader()
 {
   // Destructor
-  delete fReaderHeader;
   delete fAliHeader;
 }
 
@@ -155,9 +152,6 @@ Bool_t AliJetKineReader::FillMomentumArray(Int_t event)
       if (pdg == kNeutron || pdg == kK0Long) continue;
     } // End fast EMCAL
     
-    fMass = part->GetCalcMass();
-    fPdgC = part->GetPdgCode();
-
     // Fill momentum array
     Float_t r  = p/p0;
     Float_t px = r * part->Px(); 
index b562d62..364fbb3 100644 (file)
@@ -18,10 +18,6 @@ class AliJetKineReader : public AliJetReader
  public: 
   AliJetKineReader();
   virtual ~AliJetKineReader();
-
-  // Getters
-  Float_t GetParticleMass() const {return fMass;}        // returns mass of the Track
-  Int_t   GetParticlePdgCode() const {return fPdgC;}     // returns Pdg code
   // Setters
   Bool_t  FillMomentumArray(Int_t event);
   void    OpenInputFiles();
@@ -38,10 +34,6 @@ class AliJetKineReader : public AliJetReader
 
   AliRunLoader *fRunLoader;       //! Pointer to the run loader
   AliHeader    *fAliHeader;       //! Header
-  
-  Float_t fMass;                  //!Particle mass
-  Int_t   fPdgC;                  //!Pdg code
   ClassDef(AliJetKineReader,1)
 };
  
index ae1cc23..7f2c642 100644 (file)
 #include "AliESD.h"
 #include "AliESDtrack.h"
 
-ClassImp(AliJetMCReader)
+ClassImp(AliJetMCReader);
 
-AliJetMCReader::AliJetMCReader()
+
+AliJetMCReader::AliJetMCReader():
+    AliJetESDReader(),
+    fChainMC(0x0)
 {
   // Constructor
-  fReaderHeader = 0x0;
-  fMass = 0;
-  fPdgC = 0;
 }
 
 //____________________________________________________________________________
@@ -46,7 +46,7 @@ AliJetMCReader::AliJetMCReader()
 AliJetMCReader::~AliJetMCReader()
 {
   // Destructor
-  delete fReaderHeader;
+    delete fChainMC;
 }
 
 //____________________________________________________________________________
@@ -85,8 +85,6 @@ Bool_t AliJetMCReader::FillMomentumArray(Int_t event)
     if (pt < ptMin) continue; //check  cuts 
     p = part->P();
     e = part->Energy();
-    fMass = part->GetCalcMass();
-    fPdgC = part->GetPdgCode();
    // fill momentum array
     new ((*fMomentumArray)[goodTrack]) TLorentzVector(p.X(), p.Y(), p.Z(), e);
     goodTrack++;
index 9c754ff..dfa8f53 100644 (file)
@@ -18,7 +18,7 @@ class AliJetMCReader : public AliJetESDReader
     Bool_t FillMomentumArray(Int_t event);
     
  protected:
-    Float_t fPdgC;   // Pdg code
+    TChain* fChainMC;// Input chain
     ClassDef(AliJetMCReader,1)
 };
  
index c0a3801..fcaef90 100644 (file)
@@ -25,10 +25,10 @@ ClassImp(AliJetMCReaderHeader)
 //____________________________________________________________________________
 
 AliJetMCReaderHeader::AliJetMCReaderHeader():
- AliJetReaderHeader("AliJetMCReaderHeader") 
+    AliJetESDReaderHeader()
 {
   // Constructor
-  fDCA = 3.;
+    fDCA = 3.;
 }
 
 //____________________________________________________________________________
index d38ab4e..f971690 100644 (file)
@@ -8,9 +8,9 @@
 // Header for the MC reader in the jet analysis
 // Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
 
-#include "AliJetReaderHeader.h"
+#include "AliJetESDReaderHeader.h"
  
-class AliJetMCReaderHeader : public AliJetReaderHeader
+class AliJetMCReaderHeader : public AliJetESDReaderHeader
 {
 
  public:
index 326ca15..f31d1fc 100755 (executable)
@@ -16,7 +16,8 @@
 //------------------------------------------------------------------------
 // Jet reader base class
 // manages the reading of input for jet algorithms
-// Author: jgcn@mda.cinvestav.mx
+// Authors: jgcn@mda.cinvestav.mx
+//          magali.estienne@IReS.in2p3.fr
 //------------------------------------------------------------------------
 
 // root
@@ -32,20 +33,19 @@ ClassImp(AliJetReader)
 ////////////////////////////////////////////////////////////////////////
 
 AliJetReader::AliJetReader():
-  fChain(0),
-  fChainMC(0),
-  fMomentumArray(0),
+  fMomentumArray(new TClonesArray("TLorentzVector",2000)),
   fArrayMC(0),
-  fESD(0),
+  fFillUnitArray(new TTask("fillUnitArray","Fill unit array jet finder")),
   fReaderHeader(0),
   fSignalFlag(0),
-  fCutFlag(0)
-    
+  fCutFlag(0),
+  fUnitArray(new AliJetUnitArray[60000]),     
+  fUnitArrayNoCuts(new AliJetUnitArray[60000]),
+  fArrayInitialised(0)
 {
   // Default constructor
-  fMomentumArray = new TClonesArray("TLorentzVector",2000);
   fSignalFlag = TArrayI();
-  fCutFlag = TArrayI();
+  fCutFlag    = TArrayI();
 }
 
 ////////////////////////////////////////////////////////////////////////
@@ -53,14 +53,27 @@ AliJetReader::AliJetReader():
 AliJetReader::~AliJetReader()
 {
   // Destructor
-  delete fChain;
-  delete fChainMC;
-  delete fESD;
   if (fMomentumArray) {
       fMomentumArray->Delete();
       delete fMomentumArray;
   }
+  
+  if (fUnitArray) {
+      fUnitArray->Delete();
+      delete fUnitArray;
+  }
+  
+  if (fUnitArrayNoCuts) {
+    fUnitArrayNoCuts->Delete();
+    delete fUnitArrayNoCuts;
+  }
+
+  if (fFillUnitArray) {
+    fFillUnitArray->Delete();
+    delete fFillUnitArray;
+  }
   delete fArrayMC;
+  
 }
 
 
@@ -70,4 +83,5 @@ void AliJetReader::ClearArray()
 
 {
   if (fMomentumArray)  fMomentumArray->Clear();
+  if (fFillUnitArray)  fFillUnitArray->Clear();
 }
index 58ffc08..9bbc416 100755 (executable)
@@ -6,14 +6,24 @@
  
 // Jet reader base class
 // manages the reading of input for jet algorithms
-// Author: jgcn@mda.cinvestav.mx
-  
+// Authors: jgcn@mda.cinvestav.mx
+//          Magali Estienne <magali.estienne@IReS.in2p3.fr>  
+
 #include <TObject.h>
 #include <TChain.h>
 #include <TArrayI.h>
+#ifndef ROOT_TTask
+#include "TTask.h"
+#endif
+
+#include <AliJetUnitArray.h>
+#include <AliJetHadronCorrectionv1.h>
+
 class TTree;
+class TTask;
 class TClonesArray;
 class AliJetReaderHeader;
+class AliJetUnitArray;
 class AliESD;
 class AliJet;
 
@@ -25,15 +35,24 @@ class AliJetReader : public TObject
 
   // Getters
   virtual TClonesArray *GetMomentumArray() {return fMomentumArray;}
-  virtual Int_t GetChainEntries() {return fChain->GetEntries();} 
+
+  virtual AliJetUnitArray     *GetUnitArray() const {return fUnitArray;}  
+  virtual AliJetUnitArray     *GetUnitArrayNoCuts() const {return fUnitArrayNoCuts;}
+  
   virtual AliJetReaderHeader* GetReaderHeader() { 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 GetCutFlag(Int_t i)    const {return fCutFlag[i];}
+  virtual Int_t GetArrayInitialised() {return fArrayInitialised;}
   
   // Setters
   virtual Bool_t FillMomentumArray(Int_t) {return kTRUE;}
+  virtual void   FillUnitArrayFromTPCTracks(Int_t) {}     // temporarily not used
+  virtual void   FillUnitArrayFromEMCALHits() {}          // temporarily not used
+  virtual void   FillUnitArrayFromEMCALDigits(Int_t) {}   // temporarily not used
+  virtual void   FillUnitArrayFromEMCALClusters(Int_t) {} // temporarily not used
+  virtual void   InitUnitArray() {}
   virtual void   SetReaderHeader(AliJetReaderHeader* header) 
-    {fReaderHeader = header;}
+      {fReaderHeader = header;}
          
   // Others
   virtual void   OpenInputFiles() {}
@@ -45,17 +64,16 @@ class AliJetReader : public TObject
  protected:
   AliJetReader(const AliJetReader& rJetReader);
   AliJetReader& operator = (const AliJetReader& rhsr);
-
-  TChain                  *fChain;         // chain for reconstructed tracks
-  TChain                  *fChainMC;       // chain for mc information
-  TClonesArray            *fMomentumArray; // array of particle momenta
-  TClonesArray            *fArrayMC;       // array of mc particles
-  AliESD                  *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
-
+  TClonesArray            *fMomentumArray;    // array of particle momenta
+  TClonesArray            *fArrayMC;          // array of mc particles
+  TTask                   *fFillUnitArray;    // task list for filling the UnitArray
+  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 
+  Bool_t                   fArrayInitialised; // To check that array of units is initialised  
   ClassDef(AliJetReader,1)
 };
  
index e3733d1..796437d 100755 (executable)
@@ -27,14 +27,11 @@ ClassImp(AliJetReaderHeader)
 
 AliJetReaderHeader::AliJetReaderHeader():  
  TNamed("AliJetReaderHeader", "Jet Reader Header"),
- fNesd(0),
  fFirst(0),
  fLast(-1),
  fFiducialEtaMin(-0.9),
  fFiducialEtaMax(0.9),
  fPtCut(2.0),
- fDCA(0.),
- fTLength(0.),
  fComment("No comment"),
  fDir(""),
  fPattern("")
@@ -46,14 +43,11 @@ AliJetReaderHeader::AliJetReaderHeader():
 
 AliJetReaderHeader::AliJetReaderHeader(const char * name):  
  TNamed(name, "Jet Reader Header"),
- fNesd(0),
  fFirst(0),
  fLast(-1),
  fFiducialEtaMin(-0.9),
  fFiducialEtaMax(0.9),
  fPtCut(2.0),
- fDCA(0.),
- fTLength(0.),
  fComment("No comment"),
  fDir(""),
  fPattern("")
index 6dc9b5b..4d7c6f0 100755 (executable)
@@ -27,37 +27,32 @@ class AliJetReaderHeader : public TNamed
   virtual const char* GetPattern() {return fPattern.Data();}
   virtual Float_t     GetFiducialEtaMin() const {return fFiducialEtaMin;}
   virtual Float_t     GetFiducialEtaMax() const {return fFiducialEtaMax;}  
-  virtual Float_t GetPtCut() const {return fPtCut;}
-  Float_t GetDCA() const  {return fDCA;}       // not working so far..(always 0)
-  Float_t GetTLength() const {return fTLength;}   // not working so far.. (always 0)
-  Int_t   GetNesd() const {return fNesd;}
-  Int_t   GetNEvents() const {return fLast-fFirst;}
+  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   GetLastEvent()  const {return fLast;}
+  Int_t   GetDetector()   const {return fOption;}
+  Int_t   GetDebug()      const {return fDebug;}
   // Setters
   virtual void SetComment(const char* s) {fComment=TString(s);}
   virtual void SetPattern(const char* s) {fPattern=TString(s);}
   virtual void SetDirectory(const char* s) {fDir=TString(s);}
-  virtual void SetNumberOfESD(Int_t i=1) {fNesd = i;}
   virtual void SetFirstEvent(Int_t i=0) {fFirst=i;}
   virtual void SetLastEvent(Int_t i=-1) {fLast=i;}
   virtual void SetFiducialEta(Float_t etamin, Float_t etamax) 
       { fFiducialEtaMin = etamin; fFiducialEtaMax = etamax;}
   virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;}
-  virtual void SetDCA(Float_t dca = 0.0) {fDCA = dca;}
-  virtual void SetTLength(Float_t length = 0.0) {fTLength = length;}
 
+  virtual void SetDetector(Int_t option = 0) {fOption = option;}
+  virtual void SetDebug(Int_t debug = 0) {fDebug = debug;}
  protected:
-
-  Int_t fNesd;             // Number of ESDs to read
-  Int_t fFirst;            // First and last events analyzed
-  Int_t fLast;             // in current set of files
+  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
   Float_t fFiducialEtaMin; // Fiducial minimum eta
   Float_t fFiducialEtaMax; // Fiducial maximum eta
   Float_t fPtCut;          // pt cut
-  Float_t fDCA;            // dca cut
-  Float_t fTLength;        // track length cut
   TString fComment;        // a comment
   TString fDir;            // directory with input files
   TString fPattern;        // pattern to look for input files
index c19bb5d..3dd4cef 100755 (executable)
@@ -75,7 +75,7 @@ extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
 
 ////////////////////////////////////////////////////////////////////////
 
-void AliUA1JetFinder::FindJets()
+void AliUA1JetFinder::FindJetsTPC()
 
 {
   // initialize event, look for jets, download jet info
@@ -215,6 +215,253 @@ void AliUA1JetFinder::FindJets()
   delete percentage;
 }
 
+void AliUA1JetFinder::FindJets()
+{ // Find jets with TPC or EMCAL or TPC + EMCAL information
+  // initialize event, look for jets, download jet info
+
+  // 1) transform input to pt,eta,phi plus lego
+    AliJetUnitArray *fUnit = fReader->GetUnitArray();
+    Int_t nIn = fUnit->GetUnitEntries();
+    Int_t fOption = fReader->GetReaderHeader()->GetDetector();
+    Int_t fDebug = fReader->GetReaderHeader()->GetDebug();
+    
+    if(fDebug>1){ 
+       printf("Inside FindJets function ! \n");
+       printf("fOption : %d \n",fOption);
+    }
+    
+    if(fDebug>10){ 
+       cout << "fUnit : " << fUnit << endl;
+       printf("nIn = fUnit->GetUnitEntries() : %d \n",nIn);
+    }
+    
+    if(fDebug>10){
+       printf("     -----------------------------------------------------------------\n");
+       printf("     --- All inputs in fUnitArray ---\n");
+       for(Int_t i=0; i<nIn ; i++){
+           if(fUnit[i].GetUnitEnergy()!=0){
+               printf("     -----------------------------------------------------------------\n");
+               cout << "     |  ID : " << fUnit[i].GetUnitID() << "  |  Eta : " << fUnit[i].GetUnitEta() << "  |  Phi : " << fUnit[i].GetUnitPhi() << "  |  Energy : " << fUnit[i].GetUnitEnergy() << " |" <<endl;
+           }
+       }
+    }
+    
+    if (nIn == 0) return;
+    
+    Int_t nCandidateCut = 0;
+    Int_t nCandidateCut2 = 0;
+    Int_t nCandidate = 0;
+    for (Int_t i = 0; i < nIn; i++){
+       if(fUnit[i].GetUnitEnergy()>0. && fUnit[i].GetUnitCutFlag()==1) {
+           // Number of good tracks/digits which have passed the pt cut
+           nCandidateCut += 1;
+       }
+    if(fUnit[i].GetUnitEnergy()>0.) {
+       // Number of good tracks/digits without pt cut
+       nCandidate += 1;
+    }
+    }
+    
+    if(fDebug>=3){
+    cout << "nCandidate : " << nCandidate << endl;
+    cout << "nCandidateCut : " << nCandidateCut << endl;
+    cout << "nCandidateCut2 : " << nCandidateCut2 << endl;
+    }
+    
+  // 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;
+    
+    fJets->SetNinput(nCandidate);
+    
+    if(fDebug>3){
+       cout << "nCandidate : " << nCandidate << endl;
+       //    cout << "nMultCandidate : " << nMultCandidate << endl;
+    }
+    
+    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
+    
+    // Information extracted from fUnitArray
+    for(Int_t i=0; i<nIn; i++) 
+    {
+       if(fUnit[i].GetUnitCutFlag()==1){ 
+           etCell[i] = fUnit[i].GetUnitEnergy();
+           if (etCell[i] > 0.0) etCell[i] -= fHeader->GetMinCellEt();
+           if (etCell[i] < 0.0) etCell[i] = 0.;
+           etaCell[i] = fUnit[i].GetUnitEta();
+           phiCell[i] = fUnit[i].GetUnitPhi();
+       }
+       else {
+           etCell[i]  = 0.;
+           etaCell[i] = fUnit[i].GetUnitEta();
+           phiCell[i] = fUnit[i].GetUnitPhi();
+       }
+       
+       if(fUnit[i].GetUnitEnergy()>0.){
+           ptT[loop1]   = fUnit[i].GetUnitEnergy();
+           enT[loop1]   = fUnit[i].GetUnitEnergy();
+           etaT[loop1]  = fUnit[i].GetUnitEta();
+           phiT[loop1]  = fUnit[i].GetUnitPhi();
+           detaT[loop1] = fUnit[i].GetUnitDeta();
+           dphiT[loop1] = fUnit[i].GetUnitDphi();
+           cFlagT[loop1]= fUnit[i].GetUnitCutFlag();
+           idT[loop1]   = fUnit[i].GetUnitID();
+           loop1++;
+       }
+    }
+    
+    
+    if(fDebug > 40) // For comparison
+    {
+       for(Int_t j=0; j<nIn; j++) {
+           if(etCell[j]>0){
+               cout << "etCell[" << j << "] : " << etCell[j] << endl;
+               cout << "etaCell[" << j << "] : " << etaCell[j] << endl;
+               cout << "phiCell[" << j << "] : " << phiCell[j] << endl;
+           }
+       }
+    }
+    
+    // Run the algo. Parameters from header
+    //  Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
+    Int_t nTot      = nIn;
+    Float_t minmove = fHeader->GetMinMove();
+    Float_t maxmove = fHeader->GetMaxMove();
+    Int_t mode      = fHeader->GetMode();
+    Float_t precbg  = fHeader->GetPrecBg();
+    Int_t ierror;
+    
+    ua1_jet_finder(nIn, nTot, etCell, etaCell, phiCell, 
+                  minmove, maxmove, mode, precbg, ierror);
+    
+    // sort jets
+    Int_t * idx  = new Int_t[UA1JETS.njet];
+    TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
+    
+    if(fDebug > 20)
+    {
+       for(Int_t i = 0; i < UA1JETS.njet; i++) 
+       {
+           cout << "Number of jets found, UA1JETS.njet : " << UA1JETS.njet << endl;
+           cout << "UA1JETS.etj : " << UA1JETS.etj << endl;
+           cout << "idx[" << i << "] : " << idx[i] << endl;
+           cout << "UA1JETS.etaj[1][" << idx[i] << "] : " << UA1JETS.etaj[1][idx[i]] << endl;
+           cout << "UA1JETS.phij[1][" << idx[i] << "] : " << UA1JETS.phij[1][idx[i]] << endl;
+           cout << "UA1JETS.etj[" << idx[i] << "] : " << UA1JETS.etj[idx[i]] << endl;
+       }
+    }
+    
+    // download jet info.   
+    for(Int_t i = 0; i < UA1JETS.njet; i++) {
+       // reject events outside acceptable eta range
+       if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
+           || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
+       {
+           continue;
+       }
+       
+       Float_t px, py,pz,en,pT; // convert to 4-vector
+       px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
+       py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
+       pz = UA1JETS.etj[idx[i]] /
+           TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
+       en = TMath::Sqrt(px * px + py * py + pz * pz);
+       pT = TMath::Sqrt(px * px + py * py);
+       
+       fJets->AddJet(px, py, pz, en);
+       
+    }
+
+    // find multiplicities and relationship jet-particle
+    // find also percentage of pt from pythia
+    
+    Int_t* injet = new Int_t[nCandidate];
+    Int_t* sflag = new Int_t[nCandidate];
+    for (Int_t i = 0; i < nCandidate; i++) {injet[i]= 0;sflag[i]=0;}
+    Int_t* mult  = new Int_t[fJets->GetNJets()];
+    Int_t* ncell  = new Int_t[fJets->GetNJets()];
+    Float_t* percentage  = new Float_t[fJets->GetNJets()];
+    
+    // Instead of using etaT below, it would be interesting to use the previous fUnitArray object
+    // With the particle ID, it is possible to to have access to its physical properties and one can,
+    // for example, set if the corresponding particle is inside or outside the jet with the flag 
+    // kOutJet/kInJet, other possibilities...
+    
+    for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
+       Float_t pt_sig = 0.0;
+       mult[i] = 0;
+       ncell[i] = UA1JETS.ncellj[i];
+       for (Int_t j = 0; j < nCandidate; j++) {
+           Float_t deta = etaT[j] - fJets->GetEta(i);
+           Float_t dphi = phiT[j] - fJets->GetPhi(i);
+           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 < fHeader->GetRadius() && injet[j] == 0) {
+               injet[j] = -(i+1);
+               if(cFlagT[j] == 1 &&
+                  (etaT[j] < fHeader->GetLegoEtaMax()) &&
+                  (etaT[j] > fHeader->GetLegoEtaMin())) {
+                   injet[j] = i+1;
+                   mult[i]++;
+                   pt_sig+=enT[j];
+                   sflag[j]=1;
+               }
+           }
+           if(fDebug>10){
+               cout << "mult[" << i << "] : " << mult[i] << endl;
+               cout << "ncell[" << i << "] : " << ncell[i] << endl;
+           }
+       }
+       percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
+           ((Double_t) fJets->GetPt(i));    
+    }
+    
+    fJets->SetNCells(ncell);
+    fJets->SetPtFromSignal(percentage);
+    fJets->SetMultiplicities(mult);
+    fJets->SetInJet(injet);
+    fJets->SetEtaIn(etaT);
+    if(fDebug>10){
+       for(Int_t i=0; i<nCandidate ; i++){
+           cout << "phiT[" << i << "] : " << phiT[i] << endl;
+           cout << "etaT[" << i << "] : " << etaT[i] << endl;
+       }
+    }
+    fJets->SetPhiIn(phiT);
+    fJets->SetPtIn(enT);
+    fJets->SetEtAvg(UA1JETS.etavg);
+    delete etCell;
+    delete etaCell;
+    delete phiCell;
+    delete ncell;
+    delete cFlagT;
+    delete cClusterT;
+    delete enT;
+    delete ptT;
+    delete etaT;
+    delete phiT;
+    delete detaT;
+    delete dphiT;
+    delete injet;
+    delete idx;
+    delete mult;
+    delete percentage;
+    
+}
+
+
 ////////////////////////////////////////////////////////////////////////
 
 void AliUA1JetFinder::Reset()
index 8686b5e..afae577 100755 (executable)
@@ -28,6 +28,7 @@ class AliUA1JetFinder : public AliJetFinder
   // setters
   void SetJetHeader(AliUA1JetHeader* h) {fHeader= h;}
   // others
+  void FindJetsTPC();
   void FindJets();
   void Reset();
   void Init();
index 3d65f92..3cdee9c 100644 (file)
 #pragma link C++ class AliJetDistributions+;
 #pragma link C++ class AliUA1JetFinderV1+;
 #pragma link C++ class AliUA1JetHeaderV1+;
+#pragma link C++ class AliJetGrid+;
+#pragma link C++ class AliJetUnitArray+;
+#pragma link C++ class AliJetHadronCorrection+; 
+#pragma link C++ class AliJetHadronCorrectionv0+;
+#pragma link C++ class AliJetHadronCorrectionv1+; 
+#pragma link C++ class AliJetFillUnitArrayTracks+;
 #endif
index 47411c3..504e204 100644 (file)
@@ -7,10 +7,13 @@ SRCS =        AliJet.cxx AliJetHeader.cxx \
        AliLeading.cxx \
        AliJetProductionData.cxx AliJetProductionDataPDC2004.cxx \
        AliJetAnalysis.cxx AliJetDistributions.cxx \
-       AliUA1JetFinderV1.cxx AliUA1JetHeaderV1.cxx
+       AliUA1JetFinderV1.cxx AliUA1JetHeaderV1.cxx \
+       AliJetGrid.cxx AliJetUnitArray.cxx AliJetHadronCorrection.cxx \
+       AliJetHadronCorrectionv0.cxx AliJetHadronCorrectionv1.cxx \
+       AliJetFillUnitArrayTracks.cxx
 
 
 HDRS:= $(SRCS:.cxx=.h)
-
+EINCLUDE:= EMCAL
 DHDR= JETANLinkDef.h