Code for jet finding using TPC and EMCAL ESD information.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Nov 2006 13:45:46 +0000 (13:45 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Nov 2006 13:45:46 +0000 (13:45 +0000)
(Moved from EMCAL + updates by M. Estienne)

13 files changed:
JETAN/AliJetFillUnitArrayTracks.cxx [new file with mode: 0644]
JETAN/AliJetFillUnitArrayTracks.h [new file with mode: 0644]
JETAN/AliJetFinderTypes.h [new file with mode: 0644]
JETAN/AliJetGrid.cxx [new file with mode: 0644]
JETAN/AliJetGrid.h [new file with mode: 0644]
JETAN/AliJetHadronCorrection.cxx [new file with mode: 0644]
JETAN/AliJetHadronCorrection.h [new file with mode: 0644]
JETAN/AliJetHadronCorrectionv0.cxx [new file with mode: 0644]
JETAN/AliJetHadronCorrectionv0.h [new file with mode: 0644]
JETAN/AliJetHadronCorrectionv1.cxx [new file with mode: 0644]
JETAN/AliJetHadronCorrectionv1.h [new file with mode: 0644]
JETAN/AliJetUnitArray.cxx [new file with mode: 0644]
JETAN/AliJetUnitArray.h [new file with mode: 0644]

diff --git a/JETAN/AliJetFillUnitArrayTracks.cxx b/JETAN/AliJetFillUnitArrayTracks.cxx
new file mode 100644 (file)
index 0000000..b41ddac
--- /dev/null
@@ -0,0 +1,555 @@
+
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//======================================================================
+// *** 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
+// Author: magali.estienne@ires.in2p3.fr
+//======================================================================
+
+
+// --- Standard library ---
+#include <Riostream.h>
+
+// --- ROOT system ---
+#include <TSystem.h>
+#include <TLorentzVector.h>
+#include <TVector3.h>
+//#include "Math/Vector3D.h"
+//#include "Math/Vector3Dfwd.h"
+#include "TTask.h"
+#include <TGeoManager.h>
+#include <TMatrixD.h>
+#include <TArrayD.h>
+
+// --- AliRoot header files ---
+#include "AliJetFinder.h"
+#include "AliJetReaderHeader.h"
+#include "AliJetReader.h"
+#include "AliJetESDReader.h"
+#include "AliJetESDReaderHeader.h"
+#include "AliESD.h"
+#include "AliEMCALGeometry.h"
+#include "AliJetUnitArray.h"
+#include "AliJetFillUnitArrayTracks.h"
+#include "AliJetGrid.h"
+
+ClassImp(AliJetFillUnitArrayTracks)
+
+AliEMCALGeometry *AliJetFillUnitArrayTracks::fGeom=0; 
+
+//_____________________________________________________________________________
+AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
+  : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information")
+{
+  // constructor
+   fNIn = 0;
+  fOpt = 0;
+  fDebug = 0;
+  fNphi = 0;
+  fNeta = 0;
+  fPhiMin = 0;
+  fPhiMax = 0;
+  fEtaMin = 0;
+  fEtaMax = 0;
+  fEtaBinInTPCAcc = 0;
+  fPhiBinInTPCAcc = 0;
+  fEtaBinInEMCalAcc = 0;
+  fPhiBinInEMCalAcc = 0;
+  fNbinPhi = 0;
+}
+
+//____________________________________________________________________________
+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");
+}
+
+//____________________________________________________________________________
+void AliJetFillUnitArrayTracks::InitParameters()
+{
+  fHCorrection    = 0;     // For hadron correction
+  fHadCorr        = 0;     // For hadron correction
+  fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
+  fDebug = fReaderHeader->GetDebug();
+
+  fEtaMinCal = fGeom->GetArm1EtaMin();
+  fEtaMaxCal = fGeom->GetArm1EtaMax();
+  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, 
+                       fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
+
+  fEta   = fTPCGrid->GetArrayEta();
+  fPhi   = fTPCGrid->GetArrayPhi();
+  fIndex = fTPCGrid->GetIndexObject();
+
+  if(fDebug>20){
+    for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
+    for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
+    
+    for(Int_t i=0; i<fNphi+1; i++)
+      for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
+         (*fIndex)(i,j) << endl; }
+  } 
+  if(fDebug>1) printf("\n Parameters initiated ! \n");
+}
+
+//_____________________________________________________________________________
+AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
+{
+  // destructor
+}
+
+//_____________________________________________________________________________
+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();
+
+  if(fDebug>3){
+    cout << "nElements : " << nElements << endl;
+    cout << "nElements2 : " << nElements2 << endl;
+    cout << "fNumUnits : " << fNumUnits << endl;
+    cout << "sum : " << fNumUnits+nElements << endl;
+  }
+
+  // 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;
+           }
+
+           // 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);
+       }
+
+      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;
+           }
+
+           // Fill energy outside emcal acceptance
+           fUnitArray[fNumUnits-1+idTPC].SetUnitEnergy(unitEnergy2 + ptT[i]);
+
+           // 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)
+           {
+             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());
+                   } 
+                 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;
+                       }
+                     } 
+                 }
+                 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;
+                 }
+               }
+           }
+
+       } // end  if(fOpt==1 || option=="tpc")
+
+      delete[] enT;
+      delete[] ptT;
+      delete[] etaT;
+      delete[] phiT;
+      delete[] thetaT;
+      
+}
+
+//__________________________________________________________
+void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
+{
+  for(Int_t j=0; j<fNphi+1; j++) {
+    for(Int_t i=0; i<fNeta+1; i++) {
+
+      // TPC grid only 
+      //-------------------------------------
+      if(fGrid==0) {   
+       if(j*(fNeta+1)+i == index) {
+         eta = fEta2->At(i); 
+         phi = fPhi2->At(j);
+       }
+      }
+
+      // TPC-EMCAL grid
+      //-------------------------------------
+      Int_t ii = 0;
+      if(i==0) ii = 0;
+      if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i; 
+      if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
+
+      if(fGrid==1) {
+       if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
+         eta = fEta2->At(i);
+         phi = fPhi2->At(j);
+       }  
+
+       if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) && 
+          ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
+         if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
+         else eta = fEta2->At(i);
+         phi = fPhi2->At(j);
+       }
+
+       if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
+         eta = fEta2->At(i);
+         phi = fPhi2->At(j);
+       }
+      }
+    }
+  }
+}
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/JETAN/AliJetFillUnitArrayTracks.h b/JETAN/AliJetFillUnitArrayTracks.h
new file mode 100644 (file)
index 0000000..4dd0c6a
--- /dev/null
@@ -0,0 +1,105 @@
+#ifndef ALIJETFILLUNITARRAYTRACKS_H
+#define ALIJETFILLUNITARRAYTRACKS_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
+
+#include <TMatrixD.h>
+#include <TArrayD.h>
+
+class AliEMCALGeometry;
+class AliJetHadronCorrection;
+class AliJetReader;
+class AliJetESDReader;
+class TClonesArray;
+class AliJetUnitArray;
+//class AliJetReaderHeader;
+class AliJetReader;
+class AliJetGrid;
+
+class AliJetFillUnitArrayTracks : public TTask
+{
+ public: 
+  AliJetFillUnitArrayTracks();
+  virtual ~AliJetFillUnitArrayTracks();
+  
+  // Setter
+  void SetReaderHeader(AliJetReaderHeader *readerHeader) {fReaderHeader = readerHeader;}
+  void SetMomentumArray(TClonesArray *momentumArray) {fMomentumArray = momentumArray;}
+  void SetUnitArray(AliJetUnitArray *unitArray) {fUnitArray = unitArray;}
+  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);
+
+  // 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*);
+
+ protected:
+  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
+  AliJetHadronCorrectionv1   *fHadCorr;         // Pointer to Hadron Correction Object
+  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
+
+  AliJetReaderHeader          *fReaderHeader;   // ReaderHeader
+  TClonesArray                *fMomentumArray;  // MomentumArray
+  AliJetUnitArray             *fUnitArray;      // UnitArray
+  AliJetGrid                  *fTPCGrid;        // Define filled grid
+  AliJetGrid                  *fEMCalGrid;      // Define filled grid
+
+  // geometry info
+  static AliEMCALGeometry *fGeom;     //!
+
+  Int_t     fNphi;                    // number of points in the grid:   phi
+  Int_t     fNeta;                    //               "                 eta
+  TArrayD*  fPhi2;                    // grid points in phi
+  TArrayD*  fEta2;                    // grid points in eta
+  TArrayD*  fPhi;                     // grid points in phi
+  TArrayD*  fEta;                     // grid points in eta
+  TMatrixD* fIndex;                   // grid points in (phi,eta) 
+  TMatrixD* fParams;                  // matrix of parameters in the grid points  
+  Int_t     fGrid;                    // Select the grid acceptance you want to fill
+                                      // 0 = TPC acceptance, 1 = TPC-EMCal acceptance
+  Float_t   fPhiMin;
+  Float_t   fPhiMax;
+  Float_t   fEtaMin;
+  Float_t   fEtaMax;
+  Int_t     fEtaBinInTPCAcc;
+  Int_t     fPhiBinInTPCAcc;
+  Int_t     fEtaBinInEMCalAcc;
+  Int_t     fPhiBinInEMCalAcc;
+  Int_t     fNbinPhi;
+
+ private:
+  void SetEMCALGeometry();
+  void InitParameters();
+
+  //  Int_t fEvent;
+  
+  ClassDef(AliJetFillUnitArrayTracks,1) // Fill Unit Array with tpc and/or emcal information
+};
+
+#endif
diff --git a/JETAN/AliJetFinderTypes.h b/JETAN/AliJetFinderTypes.h
new file mode 100644 (file)
index 0000000..b83de20
--- /dev/null
@@ -0,0 +1,91 @@
+/**************************************************************************
+ * 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$ */
+
+//=========================================================================
+// Modified class for JETAN
+// Some flags have been changed
+// Some are not used for the moment
+// Author: magali.estienne@ires.in2p3.fr
+//=========================================================================
+//  Enumerated types for use in JetFinder classes
+//
+//*-- Author: Mark Horner (LBL/UCT)
+//
+
+#ifndef ALIJETFINDERBGCALCTYPE_T
+#define ALIJETFINDERBGCALCTYPE_T
+
+        typedef enum { kRatio, kCone, kConstant  
+       } AliJetFinderBGCalcType_t;
+#endif
+
+#ifndef ALIJETFINDERRESETTYPE_T
+#define ALIJETFINDERRESETTYPE_T
+
+        typedef enum {  kResetData, kResetTracks, kResetDigits, kResetParameters,
+                        kResetAll, kResetPartons, kResetParticles, kResetJets
+        } AliJetFinderResetType_t;
+#endif
+
+#ifndef  ALIJETFINDERTRACKTYPE_T
+#define  ALIJETFINDERTRACKTYPE_T
+       typedef enum {  kAllP, kEM, kCharged, kNeutral, kHadron, kChargedHadron, kNoTracks, kEMChargedPi0, kNoNeutronNeutrinoKlong
+       } AliJetFinderTrackType_t;
+#endif
+
+#ifndef  ALIJETFINDERSMEARINGTYPE_T
+#define  ALIJETFINDERSMEARINGTYPE_T
+       typedef enum {  kSmear, kEfficiency , kSmearEffic, kPerfectTracks
+       } AliJetFinderSmearingType_t;
+#endif
+       
+#ifndef  ALIJETFINDEREMCALTYPE_T
+#define  ALIJETFINDEREMCALTYPE_T
+typedef enum { kHits, kTimeCut,kNoHits, kDigits, kClusters 
+       } AliJetFinderEmcalType_t;
+#endif
+
+#ifndef  ALIJETFINDERFILETYPE_T
+#define  ALIJETFINDERFILETYPE_T
+       typedef enum {  kHijing,kPythia,kData
+       } AliJetFinderFileType_t;
+#endif
+
+#ifndef  ALIJETFINDERUA1UNITFLAGTYPE_T
+#define  ALIJETFINDERUA1UNITFLAGTYPE_T
+       typedef enum {  kInCurrentJet, kInJet, kOutJet, kBelowMinEt
+       } AliJetFinderUnitFlagType_t;
+#endif
+
+#ifndef  ALIJETFINDERUA1UNITCUTFLAGTYPE_T
+#define  ALIJETFINDERUA1UNITCUTFLAGTYPE_T
+       typedef enum {  kPtSmaller, kPtHigher
+       } AliJetFinderUnitCutFlagType_t;
+#endif
+
+#ifndef  ALIJETFINDERUA1UNITSIGNALFLAGTYPE_T
+#define  ALIJETFINDERUA1UNITSIGNALFLAGTYPE_T
+       typedef enum {  kGood, kBad
+       } AliJetFinderUnitSignalFlagType_t;
+#endif
+
+#ifndef  ALIJETFINDERUNITDETECTORFLAGTYPE_T
+#define  ALIJETFINDERUNITDETECTORFLAGTYPE_T
+       typedef enum {  kTpc, kEmcal, kAll
+       } AliJetFinderUnitDetectorFlagType_t;
+#endif
diff --git a/JETAN/AliJetGrid.cxx b/JETAN/AliJetGrid.cxx
new file mode 100644 (file)
index 0000000..f0b335a
--- /dev/null
@@ -0,0 +1,595 @@
+/**************************************************************************
+ * Copyright(c) 2001-2002, 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: AliJetGrid.cxx,v 1.0 07/09/2006 */
+
+//=========================================================================
+//  *** 07 September 2006
+//  This class allows to fill a really general grid in (eta,phi) 
+//  Two types of grid can be setted : 
+//     if fGrid==0 -> the complete acceptance of a rectangular detector acceptance is filled
+//     if fGrid==1 -> the previous grid - an other rectangular detector acceptance is filled 
+//  A (eta,phi(rad)) position can be extracted from an index value without looping over all entries
+//  An index position can be extracted from a (eta,phi(rad)) position without looping over all entries
+//
+//  How to run it :
+//  > AliJetGrid *grid = new AliJetGrid((NbinPhi-1),(NbinEta-1),PhiMin,PhiMax,EtaMin,EtaMax);
+//  > grid->SetGridType(...); // 0 or 1 for the moment
+//  > grid->SetMatrixIndexes();
+//  > grid->SetIndexIJ();
+//
+//  Author : magali.estienne@ires.in2p3.fr
+//=========================================================================
+
+// Standard headers 
+#include <Riostream.h>
+// Root headers
+#include <TMatrixD.h>
+#include <TArrayD.h>
+#include <TArrayI.h>
+// AliRoot headers
+#include "AliJetGrid.h"
+
+ClassImp(AliJetGrid)
+
+//__________________________________________________________
+AliJetGrid::AliJetGrid() {
+
+  // Default constructor
+
+  fNphi = 0;
+  fNeta = 0;
+  fPhi = 0;
+  fEta = 0;
+  fIndex = 0x0;
+  fPhiMin = 0.; // acceptance removed inside
+  fPhiMax = 0.; // acceptance removed inside
+  fEtaMin = 0.; // acceptance removed inside
+  fEtaMax = 0.; // acceptance removed inside
+  fMinPhi = 0.; // total acceptance
+  fMaxPhi = 0.; // total acceptance
+  fMinEta = 0.; // total acceptance
+  fMaxEta = 0.; // total acceptance
+  fEtaBinInTPCAcc = 0;
+  fPhiBinInTPCAcc = 0;
+  fEtaBinInEMCalAcc = 0;
+  fPhiBinInEMCalAcc = 0;
+  fNbinPhi = 0;
+
+  fDebug = 1;
+}
+
+//__________________________________________________________
+AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax) {
+
+  // Standard constructor
+
+  fNphi = nphi;
+  fNeta = neta;
+  fPhiMin = 0.; // rad - acceptance removed inside
+  fPhiMax = 0.; // rad - acceptance removed inside
+  fEtaMin = 0.; //  acceptance removed inside
+  fEtaMax = 0.; //  acceptance removed inside
+  fNbinPhi = 0;
+
+  fMaxPhi = phiMax; // total acceptance
+  fMinPhi = phiMin; // total acceptance
+  fMaxEta = etaMax; // total acceptance
+  fMinEta = etaMin; // total acceptance
+
+  fPhi  = new TArrayD(fNphi+1);
+  fEta  = new TArrayD(fNeta+1);
+  fIndexI = new TArrayI((fNeta+1)*(fNphi+1)+1);
+  fIndexJ = new TArrayI((fNeta+1)*(fNphi+1)+1);
+
+  for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = (phiMax-phiMin)/fNphi*i+phiMin;
+  for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = (etaMax-etaMin)/fNeta*i+etaMin;
+
+  if(fDebug > 3){
+    for(Int_t i=0; i<(fNphi+1); i++)  cout << (*fPhi)[i] << endl;
+    for(Int_t i=0; i<(fNeta+1); i++)  cout << (*fEta)[i] << endl;
+  }
+
+  fIndex = new TMatrixD(fNphi+1,fNeta+1);
+
+}
+
+//__________________________________________________________
+AliJetGrid::AliJetGrid(const AliJetGrid& grid):TNamed(grid) {
+
+  // Copy constructor
+
+  fNphi = grid.fNphi;
+  fNeta = grid.fNeta;
+  fPhiMin = grid.fPhiMin;
+  fPhiMax = grid.fPhiMax;
+  fEtaMin = grid.fEtaMin;
+  fEtaMax = grid.fEtaMax;
+  fEtaBinInTPCAcc = grid.fEtaBinInTPCAcc;
+  fPhiBinInTPCAcc = grid.fPhiBinInTPCAcc;
+  fEtaBinInEMCalAcc = grid.fEtaBinInEMCalAcc;
+  fPhiBinInEMCalAcc = grid.fPhiBinInEMCalAcc;
+  fNbinPhi = grid.fNbinPhi;
+  fMaxPhi = grid.fMaxPhi;
+  fMinPhi = grid.fMinPhi;
+  fMaxEta = grid.fMaxEta;
+  fMinEta = grid.fMinEta;
+
+  fPhi = new TArrayD(fNphi+1);
+  for(Int_t i=0; i<fNphi+1; i++) (*fPhi)[i] = grid.fPhi->At(i);
+  fEta = new TArrayD(fNeta+1);
+  for(Int_t i=0; i<fNeta+1; i++) (*fEta)[i] = grid.fEta->At(i);
+
+  fIndex = new TMatrixD(fNphi+1,fNeta+1);
+  for(Int_t i=0; i<fNphi+1; i++) {
+    for(Int_t j=0; j<fNeta+1; j++) (*fIndex)(i,j)=(*grid.fIndex)(i,j);
+  }
+}
+
+//__________________________________________________________
+AliJetGrid::~AliJetGrid() {
+
+  // Destructor
+  delete fPhi;
+  delete fEta;
+  delete fIndexI;
+  delete fIndexJ;
+  delete fIndex;
+}
+
+//__________________________________________________________
+void AliJetGrid::InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal) 
+{ // To set initial parameters
+
+  fPhiMin = phiMinCal; // rad
+  fPhiMax = phiMaxCal; // rad
+  fEtaMin = etaMinCal;
+  fEtaMax = etaMaxCal;
+  fNbinPhi = static_cast<int>(fPhiMin/(2*TMath::Pi()/fNphi));
+
+  // Define some binning numbers
+  if(fGrid==0){
+    for(Int_t i=0; i<fNphi+1; i++) fPhiBinInTPCAcc++;
+    fPhiBinInEMCalAcc = 0;
+    
+    for(Int_t i=0; i<fNeta+1; i++) fEtaBinInTPCAcc++;
+    fEtaBinInEMCalAcc = 0;
+  }
+
+  if(fGrid==1){
+    for(Int_t i=0; i<fNphi+1; i++) {
+      fPhiBinInTPCAcc++;
+      if(fPhi->At(i) >= fPhiMin &&
+        fPhi->At(i) <= fPhiMax)
+       fPhiBinInEMCalAcc++;
+    }
+    for(Int_t i=0; i<fNeta+1; i++) {
+      fEtaBinInTPCAcc++;
+      if((fEta->At(i) >= fEtaMin) &&
+        (fEta->At(i) <= fEtaMax))
+       fEtaBinInEMCalAcc++;
+    }
+  }
+
+}
+
+//__________________________________________________________
+TArrayD* AliJetGrid::GetArrayEta() 
+{ // Returns an array with the eta points
+
+  return fEta;
+
+}
+
+//__________________________________________________________
+TArrayD* AliJetGrid::GetArrayPhi() 
+{ // Returns an array with the phi points
+
+  return fPhi;
+
+}
+
+//__________________________________________________________
+TMatrixD* AliJetGrid::GetIndexObject()
+{ // Returns a pointer to the matrix
+
+  return fIndex;
+
+}
+
+//__________________________________________________________
+void AliJetGrid::GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, Float_t &maxphi, 
+                               Float_t &mineta, Float_t &maxeta)
+{ // Returns EMCAL acceptance initially setted
+
+  nphi = fNphi;
+  neta = fNeta;
+  minphi = fPhiMin;
+  maxphi = fPhiMax;
+  mineta = fEtaMin;
+  maxeta = fEtaMax;
+}
+
+//__________________________________________________________
+void AliJetGrid::GetBinParam(Int_t &phibintpc, Int_t &etabintpc, 
+                               Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi)
+{ // Returns number of bins in TPC and EMCAL geometry
+
+  etabintpc = fEtaBinInTPCAcc;
+  phibintpc = fPhiBinInTPCAcc;
+  etabinemc = fEtaBinInEMCalAcc;
+  phibinemc = fPhiBinInEMCalAcc;
+  nbinphi = fNbinPhi;
+}
+
+//__________________________________________________________
+Int_t AliJetGrid::GetIndexFromEtaPhi(Double_t phi,Double_t eta) const 
+{ // Tells the index value of a corresponding (eta,phi) real position
+  // Loop over all entries -> takes time. 
+  // Used one time at the begining to fill the grids
+
+  /*   this is how bins are numbered
+   
+       in all TPC        or      in TPC - EMCAL
+                                ... ... ... ..
+                                ---------------
+          ...  ... .            10 |   |   | 11
+          ---+---+---           ---------------
+    ^      6 | 7 | 8     or      8 |   |   | 9 
+    |     ---+---+---           --------------- 
+           3 | 4 | 5             4 | 5 | 6 | 7
+   phi    ---+---+---           ---------------
+           0 | 1 | 2             0 | 1 | 2 | 3  
+
+       
+             eta ->
+  */
+
+  Int_t etaBin=0,phiBin=0,absID=0;
+  Int_t etaBin2=0,etaBin3=0;
+
+  // Fill all the grid in eta/phi (all TPC acceptance)
+  //-----------------------------------------------------
+  if(fGrid==0){ 
+    if(eta <= fEta->At(0)) {
+      etaBin = 0;
+    } else if(eta >= fEta->At(fNeta)) {
+      etaBin = fNeta;
+    } else {
+      for(Int_t i=0; i<fNeta+1; i++) {
+       if(eta < fEta->At(i)) {
+         etaBin = i-1;
+         break;
+       } 
+      }
+    }
+    if(phi <= fPhi->At(0)) {
+      phiBin = 0;
+    } else if(phi >= fPhi->At(fNphi)) {
+      phiBin = fNphi;
+    } else {
+      for(Int_t i=0; i<fNphi+1; i++) {
+       if(phi < fPhi->At(i)) {
+         phiBin = i-1;
+         break;
+       } 
+      }
+    }
+    
+    // Calculate absolute id
+    absID = phiBin*(fNeta+1) + etaBin;
+
+  }
+
+  // Fill the grid but do not count id in EMCal acceptance
+  //------------------------------------------------------
+  if((eta >= fEtaMin && eta <= fEtaMax) &&
+     (phi >= fPhiMin && phi <= fPhiMax)){
+    etaBin = etaBin2 = etaBin3 = -1;
+    phiBin = -1;
+  }  
+
+  if(fGrid == 1){ 
+    if(phi<fPhiMin) {
+      if(eta <= fEta->At(0)) {
+       etaBin = 0;
+      } else if(eta >= fEta->At(fNeta)) {
+       etaBin = fNeta;
+      } else {
+       for(Int_t i=0; i<fNeta+1; i++) {
+         if(eta < fEta->At(i)) {
+           etaBin = i-1;
+           break;
+         } 
+       }
+      }
+    }
+    else {
+      if((phi>=fPhiMin) && (phi<=fPhiMax)) {
+       if(eta <= fEta->At(0)) {
+         etaBin2 = 0;
+       } else if(eta >= fEta->At(fNeta)) {
+         etaBin2 = fNeta-fEtaBinInEMCalAcc;
+       } else {
+         for(Int_t i=0; i<fNeta+1; i++) {
+           if(eta < fEta->At(i) && ((fEta->At(0)<eta && eta<fEtaMin) || 
+                                    (fEtaMax<eta && eta<fEta->At(fNeta)))) {
+             //                 cout << "i : " << i << endl;
+             if(eta<fEtaMin)
+               etaBin2 = i-1;
+             else etaBin2 = i-1-fEtaBinInEMCalAcc;
+             break;
+           } 
+         }
+       }
+      }
+      else {
+       if(eta <= fEta->At(0)) {
+         etaBin3 = 0;
+       } else if(eta >= fEta->At(fNeta)) {
+         etaBin3 = fNeta;
+       } else {
+         for(Int_t i=0; i<fNeta+1; i++) {
+           if(eta < fEta->At(i)) {
+             etaBin3 = i-1;
+             break;
+           } 
+         }
+       }
+      }
+    }
+    
+    if(phi <= fPhi->At(0)) {
+      phiBin = 0;
+    } else if(phi >= fPhi->At(fNphi)) {
+      phiBin = fNphi;
+    } else {
+      for(Int_t i=0; i<fNphi+1; i++) {
+       if(phi < fPhi->At(i)) {
+         phiBin = i-1;
+         break;
+       } 
+      }
+    }
+
+    // Calculate absID
+    //-------------------------------------------------------
+    if(phi<fPhiMin)     
+      absID = phiBin*(fNeta+1) + etaBin;
+    if(phi>=fPhiMin && phi<=fPhiMax) {
+      if(eta >= fEtaMin && eta <= fEtaMax) absID = -1;
+      else{
+       absID = (fNbinPhi+1)*(fNeta+1)
+         + (phiBin-fNbinPhi-1)*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
+         + etaBin2;
+      }
+    }
+    if(phi>fPhiMax)
+      absID = (fNbinPhi+1)*(fNeta+1)
+       + fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc)) 
+       + (phiBin-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)
+       + etaBin3;
+    
+  } // END OPTION==1    
+  
+  return absID;
+  
+}
+
+//__________________________________________________________
+void AliJetGrid::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
+{ // Get (eta,phi) position for a given index BUT loop over all entries (takes time)
+
+  for(Int_t j=0; j<fNphi+1; j++) {
+    for(Int_t i=0; i<fNeta+1; i++) {
+
+      // TPC grid only 
+      //-------------------------------------
+      if(fGrid==0) {   
+       if(j*(fNeta+1)+i == index) {
+         eta = fEta->At(i); 
+         phi = fPhi->At(j);
+       }
+      }
+
+      // TPC-EMCAL grid
+      //-------------------------------------
+      Int_t ii = 0;
+      if(i==0) ii = 0;
+      if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i; 
+      if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
+
+      if(fGrid==1) {
+       if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
+         eta = fEta->At(i);
+         phi = fPhi->At(j);
+       }  
+
+       if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) && 
+          ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
+         if(ii==0) {Int_t ind = 0; eta = fEta->At(ind);}
+         else eta = fEta->At(i);
+         phi = fPhi->At(j);
+       }
+
+       if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && 
+          ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))
+           +(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
+         eta = fEta->At(i);
+         phi = fPhi->At(j);
+       }
+      }
+    }
+  }
+}
+
+//__________________________________________________________
+Int_t AliJetGrid::GetIndex(Double_t phi, Double_t eta) 
+{ // Get index value for a (eta,phi) position - Direct value
+
+  Int_t ieta = GetIndexJFromEta(eta);
+  Int_t iphi = GetIndexIFromPhi(phi);
+
+  Int_t index = GetMatrixIndex(iphi,ieta);
+
+  if(fDebug>10){
+    cout << "(phi,eta) : " << phi << ", " << eta << endl;
+    cout << "index : " << index << endl;
+  }
+  return index;
+
+}
+
+//__________________________________________________________
+Int_t AliJetGrid::GetIndexJFromEta(Double_t eta)
+{ // Get eta id 
+  // Eta discretized
+
+  Int_t idEta =0;
+  Double_t temp = (eta+fMaxEta)/(fMaxEta-fMinEta)*fNeta;
+  if(fDebug>20)
+    {
+      cout << "eta : " << eta << endl;
+      cout << "fMaxEta : " << fMaxEta << endl;
+      cout << "fMinEta : " << fMinEta << endl;
+      cout << "fNeta : " << fNeta << endl;
+      cout << "temp eta before cast : " << temp << endl;
+    }
+  idEta = static_cast<Int_t>(temp+0.5);
+  if(fDebug>20) cout << "temp eta after cast : " << idEta << endl;
+  return idEta;
+}
+//__________________________________________________________
+Int_t AliJetGrid::GetIndexIFromPhi(Double_t phi)
+{ // Get phi id
+  // Phi discretized
+
+  Int_t idPhi = 0;
+  Double_t temp = 0.;
+  if(fMinPhi==0) temp = phi/(fMaxPhi-fMinPhi)*fNphi;
+  else temp = (phi-fMinPhi)/(fMaxPhi-fMinPhi)*fNphi;
+
+  if(fDebug>20)
+    {
+      cout << "phi : " << phi << endl;
+      cout << "fMaxPhi : " << fMaxPhi << endl;
+      cout << "fMinPhi : " << fMinPhi << endl;
+      cout << "fNphi : " << fNphi << endl;
+      cout << "temp phi before cast : " << temp << endl;
+    }
+  idPhi = static_cast<Int_t>(temp+0.5);
+  if(fDebug>20) cout << "temp phi after cast : " << idPhi << endl;
+  return idPhi;
+
+
+}
+
+//__________________________________________________________
+void AliJetGrid::SetMatrixIndex(Int_t i,Double_t par) 
+{ // Allows to set parameters using only one index (if fGrid==0) !!
+  // Not used !
+
+  Int_t iphi = (Int_t)i/fNeta;
+  Int_t ieta = i-iphi*fNeta;
+  SetMatrixIndex(iphi,ieta,par);
+
+  return;
+}
+
+//__________________________________________________________
+void AliJetGrid::SetMatrixIndexes() 
+{ // Fill the final matrix object with the corresponding index in eta/phi
+
+  for(Int_t i=0; i<fNphi+1; i++){
+    for(Int_t j=0; j<fNeta+1; j++){
+      (*fIndex)(i,j) = GetIndexFromEtaPhi(fPhi->At(i),fEta->At(j))+1;
+      if(fDebug>2){
+       cout << "(*fIndex)(" << i << "," << j << ") : " << (*fIndex)(i,j) << 
+         ", phi : " << fPhi->At(i) << ", eta : " << fEta->At(j) << endl;
+      }
+    }
+  }
+  printf("##########################################################\n");
+  printf("TMatrix object filled !\n");  
+  printf("Size of the object : phi x eta = (fNphi+1) x (fNeta+1) = %d\n",(fNphi+1)*(fNeta+1));
+  printf("##########################################################\n");
+}
+
+//__________________________________________________________
+void AliJetGrid::SetIndexIJ()
+{ // 
+
+  for(Int_t i=0; i<fNphi+1; i++){
+    for(Int_t j=0; j<fNeta+1; j++){
+      Int_t id = static_cast<Int_t>((*fIndex)(i,j));
+
+      if(id!=-1)
+       {
+         (*fIndexI)[id] = i;
+         (*fIndexJ)[id] = j;
+       }
+    }
+  }
+
+  printf("##########################################################\n");
+  printf("     In SetIndexIJ - Grid indexes setted !\n");
+  printf("##########################################################\n");
+}
+
+//__________________________________________________________
+void AliJetGrid::GetIJFromIndex(Int_t index, Int_t i, Int_t j)
+{ // returns i position id of eta and j position id of phi for a given grid index
+  i = (*fIndexI)[index];
+  j = (*fIndexJ)[index];
+}
+
+//__________________________________________________________
+void AliJetGrid::GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta)
+{ // returns eta, phi values for a given grid index
+
+  phi = fPhi->At((*fIndexI)[index]);
+  eta = fEta->At((*fIndexJ)[index]);
+}
+
+//__________________________________________________________
+Int_t AliJetGrid::GetNEntries()
+{ // Returns the number of entries of the grid
+
+  if(fDebug>20){
+    cout << "fMaxPhi : " << fMaxPhi << endl;
+    cout << "fMaxEta : " << fMaxEta << endl;
+  }
+
+  Int_t indexNum = GetIndex(fMaxPhi,fMaxEta);
+  if(fDebug>20) cout << "indexNum : " << indexNum << endl;
+  return indexNum;
+
+}
+
+//__________________________________________________________
+Int_t AliJetGrid::GetNEntries2()
+{ // Returns the number of entries of the grid
+
+  Int_t indexNum = GetIndex(fMaxPhi-1.,fMaxEta-0.5);
+    if(fDebug>20) cout << "indexNum : " << indexNum << endl;
+  return indexNum;
+
+}
+
+
+
+
+
+
diff --git a/JETAN/AliJetGrid.h b/JETAN/AliJetGrid.h
new file mode 100644 (file)
index 0000000..f522166
--- /dev/null
@@ -0,0 +1,107 @@
+#ifndef ALIJETGRID_H
+#define ALIJETGRID_H
+/* Copyright(c) 2001-2002, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// Class description :
+//
+// Author : Magali Estienne, IPHC Strasbourg - e-mail: magali.estienne@ires.in2p3.fr
+//
+// --- Standard library ---
+#include <Riostream.h>
+//-- Root headers ---
+#include <TNamed.h>
+#include <TMatrixD.h>
+#include <TArrayD.h>
+#include <TArrayI.h>
+//-------------------
+
+class AliJetGrid : public TNamed {
+
+ public:
+  
+  AliJetGrid();
+  AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t etaMin,Double_t phiMax,Double_t etaMax);
+  AliJetGrid(const AliJetGrid& grid);
+  virtual ~AliJetGrid();
+
+  // Getter
+  TArrayD*  GetArrayEta();
+  TArrayD*  GetArrayPhi();
+  TMatrixD* GetIndexObject();
+  void      InitParams(Double_t phiMinCal,Double_t phiMaxCal,Double_t etaMinCal,Double_t etaMaxCal);
+  Int_t     GetIndexFromEtaPhi(Double_t phi,Double_t eta) const;
+  void      GetEtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi);
+  Int_t     GetBinsEta() const {return fNeta+1;}
+  Int_t     GetBinsPhi() const {return fNphi+1;}
+  Int_t     GetTotBins() const {return (fNeta+1)*(fNphi+1);}
+  Int_t     GetPointsEta() const {return fNeta;}
+  Int_t     GetPointsPhi() const {return fNphi;}
+  Int_t     GetTotPoints() const {return fNeta*fNphi;}
+  Int_t     GetMatrixIndex(Int_t iphi,Int_t ieta) {return static_cast<Int_t>((*fIndex)(iphi,ieta));}
+  Int_t     GetIndexIFromPhi(Double_t phi);
+  Int_t     GetIndexJFromEta(Double_t eta);
+  Int_t     GetIndex(Double_t phi,Double_t eta);
+  void      GetAccParam(Int_t &nphi, Int_t &neta, Float_t &minphi, 
+                       Float_t &maxphi, Float_t &mineta, Float_t &maxeta);
+  void      GetBinParam(Int_t &phibintpc, Int_t &etabintpc, 
+                       Int_t &phibinemc, Int_t &etabinemc, Int_t &nbinphi);
+  void      GetIJFromIndex(Int_t index, Int_t i, Int_t j);
+  void      GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta);
+  Int_t     GetNEntries();
+  Int_t     GetNEntries2();
+  Int_t     GetDeta() {return static_cast<Int_t>((fEtaMax-fEtaMin)/fNeta); 
+    if(fDebug>21) cout << "static_cast<Int_t>((fEtaMax-fEtaMin)/fNeta) : " << 
+      static_cast<Int_t>((fEtaMax-fEtaMin)/fNeta);}
+  Int_t     GetDphi() {return static_cast<Int_t>((fPhiMax-fPhiMin)/fNphi); 
+    if(fDebug>21) cout << "static_cast<Int_t>((fPhiMax-fPhiMin)/fNphi) : " << 
+      static_cast<Int_t>((fPhiMax-fPhiMin)/fNphi);}
+  Int_t     GetGridType() {return fGrid;}
+
+  // Setter
+  void      SetEtaRange(Double_t etaMin, Double_t etaMax) {fEtaMin = etaMin; fEtaMax = etaMax;}
+  void      SetPhiRange(Double_t phiMin, Double_t phiMax) {fPhiMin = phiMin; fPhiMax = phiMax;}
+  void      SetNeta(Int_t neta) {fNeta = neta;}
+  void      SetNphi(Int_t nphi) {fNphi = nphi;}
+  void      SetMatrixIndexes();
+  void      SetMatrixIndex(Int_t i,Double_t par);
+  void      SetMatrixIndex(Int_t iphi,Int_t ieta,Double_t par) {
+                                 (*fIndex)(iphi,ieta)=par; return; }
+  void      SetGridType(Int_t type) {fGrid = type;}
+  void      SetIndexIJ();
+
+ private:
+  Int_t     fGrid;              // Close the type of grid you want to fill
+                                // 0 = grid in tpc acceptance, 1 = grid in (tpc-emcal) acceptance  
+  Int_t     fNphi;              // number of points in the grid in phi
+  Int_t     fNeta;              //               "                 eta
+  TArrayD*  fPhi;               // grid points in phi
+  TArrayD*  fEta;               // grid points in eta
+  TMatrixD* fIndex;             // matrix of indexes in the grid points  
+  TArrayI*  fIndexI;            // grid points in phi
+  TArrayI*  fIndexJ;            // grid points in eta
+  Double_t  fPhiMin;            // phi acceptance min
+  Double_t  fPhiMax;            // phi acceptance max
+  Double_t  fEtaMin;            // eta acceptance min 
+  Double_t  fEtaMax;            // eta acceptance max
+  Int_t     fEtaBinInTPCAcc;    // number of points in TPC acceptance in eta
+  Int_t     fPhiBinInTPCAcc;    // number of points in TPC acceptance in phi
+  Int_t     fEtaBinInEMCalAcc;  // number of points in EMCal acceptance in eta
+  Int_t     fPhiBinInEMCalAcc;  // number of points in EMCal acceptance in phi
+  Int_t     fNbinEta;
+  Int_t     fNbinPhi;
+  Double_t  fMaxPhi;
+  Double_t  fMinPhi;
+  Double_t  fMaxEta;
+  Double_t  fMinEta;
+  Int_t     fDebug;
+
+  ClassDef(AliJetGrid,1) // Parameters used by AliTPCtrackerParam 
+};
+
+
+#endif
+
+
+
+
diff --git a/JETAN/AliJetHadronCorrection.cxx b/JETAN/AliJetHadronCorrection.cxx
new file mode 100644 (file)
index 0000000..efd5a20
--- /dev/null
@@ -0,0 +1,35 @@
+/**************************************************************************
+ * Copyright(c) 1998-2002, 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$ */
+
+//===============================================================
+// To be modified for hadron correction using particle ID and track merging
+// Author : magali.estienne@ires.in2p3.fr
+//===============================================================
+//--
+//--
+//--  Hadron correction base class
+//--
+//--
+//--
+
+// --- AliRoot header files ---
+#include "AliJetHadronCorrection.h"
+
+ClassImp(AliJetHadronCorrection)
+
+AliJetHadronCorrection::AliJetHadronCorrection(const char *name,const char *title) 
+:TNamed(name,title) { }
diff --git a/JETAN/AliJetHadronCorrection.h b/JETAN/AliJetHadronCorrection.h
new file mode 100644 (file)
index 0000000..3a7635b
--- /dev/null
@@ -0,0 +1,27 @@
+#ifndef ALIJETHADRONCORRECTION_H
+#define ALIJETHADRONCORRECTION_H
+/* Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+//                  
+//*-- Author: Aleksei Pavlinov (WSU)
+// This pure abstract class which defines only interface
+// 
+#include "TNamed.h"
+
+class AliJetHadronCorrection : public TNamed {
+
+  public:
+  AliJetHadronCorrection(const char *name="name", const char *title="title");
+  virtual ~AliJetHadronCorrection() {} 
+
+  // Add for particle
+  virtual Double_t GetEnergy(Double_t pmom, Double_t eta, Int_t gid)=0; 
+
+  ClassDef(AliJetHadronCorrection,1) // Hadron correction for EMC (abstract class)
+};
+
+#endif // ALIJETHADRONCORRECTION_H
diff --git a/JETAN/AliJetHadronCorrectionv0.cxx b/JETAN/AliJetHadronCorrectionv0.cxx
new file mode 100644 (file)
index 0000000..6fe34e8
--- /dev/null
@@ -0,0 +1,85 @@
+/**************************************************************************
+ * Copyright(c) 1998-2002, 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$ */
+
+//===============================================================
+// To be modified for hadron correction using particle ID and track merging
+// Author : magali.estienne@ires.in2p3.fr
+//===============================================================
+
+
+// --- AliRoot header files ---
+#include "AliJetHadronCorrectionv0.h"
+
+const Int_t maxVariant = 8;  // size eta grid
+const Int_t nVec = 10;       // size momentum grid
+const Int_t nPol = 4;        // number coefficients of polinom
+static Double_t etaGrid[maxVariant]={ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.67};
+static Double_t pGrid[nVec]={ 0.2, 0.5, 1.0, 2.0, 3.0, 5.0,10.0,15.0,25.0,40.0};
+// c[][] - first index for eta, second for momentum 
+static Double_t c[maxVariant][nPol] ={
+    {1.305705e-01, 3.725653e-01, -1.219962e-02, 1.806235e-04},
+    {1.296153e-01, 3.827408e-01, -1.238640e-02, 1.822804e-04},
+    {1.337690e-01, 3.797454e-01, -1.245227e-02, 1.848243e-04},
+    {1.395796e-01, 3.623994e-01, -9.196803e-03, 1.243278e-04},
+    {1.457184e-01, 3.753655e-01, -1.035324e-02, 1.473447e-04},
+    {1.329164e-01, 4.219044e-01, -1.310515e-02, 1.948883e-04},
+    {8.136581e-02, 4.646087e-01, -1.531917e-02, 2.274749e-04},
+    {1.119836e-01, 4.262497e-01, -1.160125e-02, 1.628738e-04} };
+
+ClassImp(AliJetHadronCorrectionv0)
+
+AliJetHadronCorrectionv0* AliJetHadronCorrectionv0::fHadrCorr = 0;
+
+AliJetHadronCorrectionv0::AliJetHadronCorrectionv0(const char *name,const char *title) 
+                           :AliJetHadronCorrection(name, title)
+{
+  fHadrCorr = this;
+}
+
+AliJetHadronCorrectionv0*
+AliJetHadronCorrectionv0::Instance()
+{
+  fHadrCorr = new AliJetHadronCorrectionv0();
+  return fHadrCorr;
+}
+
+Double_t 
+AliJetHadronCorrectionv0::GetEnergy(Double_t pmom, Double_t eta, Int_t /*gid*/)
+{
+  Int_t iEta=0; // index 
+  Double_t etaw = TMath::Abs(eta);
+  if(etaw > etaGrid[maxVariant-1]) etaw = etaGrid[maxVariant-1];
+  for(Int_t i=0; i<maxVariant; i++) if(eta>=etaGrid[i]) {iEta = i; break;}
+
+  Double_t e[2], y, pw = pmom;
+  if(pmom > pGrid[nVec-1]) pw = pGrid[nVec-1];
+  for(Int_t i=0; i<2; i++){ // e for two eta value
+    e[i] = c[iEta][0];
+    y = 1.;
+    for(Int_t j=1; j<nPol; j++){
+      y *= pw;
+      e[i] += c[iEta][j]*y;
+    }
+    if(i==0) iEta ++;
+  }
+
+  Double_t deta = etaGrid[iEta] - etaGrid[iEta-1];
+  Double_t a = (e[1] - e[0])/deta; // slope 
+  Double_t energy = e[0] + a*(eta-etaGrid[iEta-1]);
+  return energy;
+}
diff --git a/JETAN/AliJetHadronCorrectionv0.h b/JETAN/AliJetHadronCorrectionv0.h
new file mode 100644 (file)
index 0000000..68d53eb
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef ALIJETHADRONCORRECTIONV0_H
+#define ALIJETHADRONCORRECTIONV0_H
+/* Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+//                  
+//*-- Author: Aleksei Pavlinov (WSU)
+#include "AliJetHadronCorrection.h"
+
+class AliJetHadronCorrectionv0: public AliJetHadronCorrection {
+
+  private:
+  static AliJetHadronCorrectionv0* fHadrCorr;
+
+  protected:
+  AliJetHadronCorrectionv0(const char *name="HadronCorrectionv0", const char *title="title");
+
+  public:
+  static  AliJetHadronCorrectionv0* Instance();
+  virtual Double_t GetEnergy(Double_t pmom, Double_t eta, Int_t gid); 
+  Double_t GetEnergy(Double_t pmom, Double_t eta) 
+  {return GetEnergy(pmom,eta,7);}
+
+  virtual ~AliJetHadronCorrectionv0() {}
+
+  ClassDef(AliJetHadronCorrectionv0,1) // Hadron correction for EMC (version for MDC)
+};
+
+#endif // ALIJETHADRONCORRECTIONV0_H
diff --git a/JETAN/AliJetHadronCorrectionv1.cxx b/JETAN/AliJetHadronCorrectionv1.cxx
new file mode 100644 (file)
index 0000000..939c81f
--- /dev/null
@@ -0,0 +1,145 @@
+/**************************************************************************
+ * Copyright(c) 1998-2002, 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.                  *
+ **************************************************************************/
+
+//===============================================================
+// To be modified for hadron correction using particle ID and track merging
+// Author : magali.estienne@ires.in2p3.fr
+//===============================================================
+// Author : Mark Horner (LBL/UCT)
+
+// --- Standard library ---
+#include "Riostream.h"
+#include "TMath.h"
+
+// --- AliRoot header files ---
+//#include "../EMCAL/AliJetGeometry.h"
+#include "AliEMCALGeometry.h"
+#include "AliJetHadronCorrectionv1.h"
+
+ClassImp(AliJetHadronCorrectionv1)
+
+Double_t AliJetHadronCorrectionv1::fgParLookup[HCPARAMETERS][HCPARAMETERSETS] = 
+{  
+    {-2.82271e-4 , -2.39954e-4},   
+    {2.50796e-2  , 2.07172e-2},   
+    {1.02861e-3  , 1.48576e-3},   
+    {2.11539e-2  , -1.38473-2},   
+    {2.27003e-2  , 2.78252e-2},   
+    {1.65078e-6  , 1.51821e-6}
+};
+
+AliJetHadronCorrectionv1* AliJetHadronCorrectionv1::fgHadrCorr = 0;
+
+void AliJetHadronCorrectionv1::SetGeometry(AliEMCALGeometry *geometry)
+{
+  // Initialise EMCAL geometry
+    if (!geometry)
+    {
+       SetParameters();
+       fSamplingFraction = geometry->GetSampling();
+       cout<<"Setting the sampling fraction to :"<<fSamplingFraction<<endl; 
+    }else
+    {
+        SetParameters(geometry->GetName());
+       fSamplingFraction = geometry->GetSampling();
+       cout<<"Setting the sampling fraction to :"<<fSamplingFraction<<endl; 
+    }
+    return;    
+}      
+
+/*
+###########################################################################
+###########################################################################
+
+This will have to be modified with the inclusion of the new EMCAL geometry
+That means I guess a new study of the hadron correction...
+
+Geometry : SHISH...
+
+###########################################################################
+###########################################################################
+*/
+
+
+
+       
+void AliJetHadronCorrectionv1::SetGeometry(TString name,Double_t fs)
+{
+  // Initialise EMCAL geometry
+  cout << "Setting sampling fraction to "<<fSamplingFraction<<endl;    
+   fSamplingFraction = fs;     
+  if ( name == ""              ||
+       name == "EMCAL_5655_21" ||
+       name == "EMCALArch1a"   ||
+       name == "EMCALArch1aN"  ||
+       name == "G56_2_55_19"   ||
+       name == "G56_2_55_19_104_14" )
+  { // set parameters to this hadron correction
+     cout<<"HC parameters!"<<endl;
+     for (Int_t i=0;i<6;i++)
+     {
+          fPar[i] = fgParLookup[i][0];  
+          cout <<fPar[i]<<endl;
+     }
+  }else if( name == "EMCAL_6564_21" ||  
+           name == "G65_2_64_19" )
+  {      
+      cout<<"HC parameters!"<<endl;
+     for (Int_t i=0;i<6;i++)
+     {
+          fPar[i] = fgParLookup[i][1];  
+          cout <<fPar[i]<<endl;
+     }
+  }else
+  {
+    printf("Geometry not defined in hadron correction\n"); 
+  }      
+       
+}      
+
+       
+AliJetHadronCorrectionv1::AliJetHadronCorrectionv1(const char *name,const char *title) 
+                           :AliJetHadronCorrection(name, title)
+{
+  fgHadrCorr = this;
+}
+
+/*
+AliJetHadronCorrectionv1::AliJetHadronCorrectionv1(const char *name,const char *title,AliJetGeometry *geometry)
+{
+
+  fHadrCorr = this;
+  SetGeometry(geometry);  
+       
+}      
+*/
+
+AliJetHadronCorrectionv1*
+AliJetHadronCorrectionv1::Instance()
+{
+  // return pointer to global instance. Instantiate if needed
+  if (! fgHadrCorr) new AliJetHadronCorrectionv1();
+  return fgHadrCorr;
+}
+
+Double_t 
+AliJetHadronCorrectionv1::GetEnergy(Double_t pmom, Double_t eta, Int_t /*gid*/)
+{
+  // Return parametrised energy response
+  Double_t etai = TMath::Abs(eta); 
+  Double_t value =  fPar[5]*pmom*pmom*pmom+ fPar[0]*pmom*pmom+fPar[1]*pmom +fPar[2]*pmom*etai +fPar[3]*etai + fPar[4];
+  return fSamplingFraction*value;
+   
+}
diff --git a/JETAN/AliJetHadronCorrectionv1.h b/JETAN/AliJetHadronCorrectionv1.h
new file mode 100644 (file)
index 0000000..6da5fd2
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALIJETHADRONCORRECTIONV1_H
+#define ALIJETHADRONCORRECTIONV1_H
+/* Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+
+//                  
+//*-- Author: Mark Horner (LBL/UCT)
+//
+#include "AliJetHadronCorrection.h"
+
+#define HCPARAMETERS    6
+#define HCPARAMETERSETS 2
+
+class AliEMCALGeometry;
+
+class AliJetHadronCorrectionv1 : public AliJetHadronCorrection 
+{
+ public:
+  static  AliJetHadronCorrectionv1* Instance();
+  virtual Double_t GetEnergy(Double_t pmom, Double_t eta, Int_t gid); 
+  Double_t GetEnergy(Double_t pmom, Double_t eta) {return GetEnergy(pmom,eta,7);}
+  
+  void SetGeometry(TString name, Double_t fs = 1.); 
+  virtual ~AliJetHadronCorrectionv1() {}
+
+ protected:
+  AliJetHadronCorrectionv1(const char *name="HadronCorrectionv1", const char *title="Hadron Correction");
+  //    AliJetHadronCorrectionv1(const char *name="HadronCorrectionv1", const char *title="Hadron Correction",AliJetGeometry *geometry = NULL);
+  void SetGeometry(AliEMCALGeometry *geometry);
+    
+ private:
+  void SetParameters(TString name = "") {Warning("SetParameter","Dummy method with argument %s",name.Data());}
+    
+  static AliJetHadronCorrectionv1* fgHadrCorr;  // Pointer to global instance (singleton)
+  static Double_t fgParLookup[HCPARAMETERS][HCPARAMETERSETS]; // Global array with parameters for hadronic response
+  Double_t fPar[6];
+  Float_t  fSamplingFraction;  // Sampling fraction
+      
+  ClassDef(AliJetHadronCorrectionv1,2) // Hadron correction for EMC (version for MDC)
+};
+
+#endif // ALIJETHADRONCORRECTIONV1_H
diff --git a/JETAN/AliJetUnitArray.cxx b/JETAN/AliJetUnitArray.cxx
new file mode 100644 (file)
index 0000000..2150d79
--- /dev/null
@@ -0,0 +1,72 @@
+/**************************************************************************
+ * 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$ */
+
+//_________________________________________________________________________
+//  Unit used by UA1 algorithm
+// --
+//*-- Author: Sarah Blyth (LBL/UCT)
+// --
+// Revised Version for JETAN - 30/03/2006
+// -- Magali Estienne (IReS)
+
+#include "AliJetUnitArray.h"
+
+
+ClassImp(AliJetUnitArray)
+
+AliJetUnitArray::AliJetUnitArray()
+{
+  // Default constructor
+  fUnitEnergy         = 0.0;
+  fUnitEta            = 0.0;
+  fUnitPhi            = 0.0;
+  fUnitID             = 0;
+  fUnitFlag           = kOutJet;
+}  
+
+AliJetUnitArray::~AliJetUnitArray()
+{
+  // Destructor 
+}
+       
+Bool_t AliJetUnitArray::operator>(AliJetUnitArray unit) const
+{
+  // Greater than operator used by sort
+  if( fUnitEnergy > unit.GetUnitEnergy())
+    return kTRUE;
+  else 
+    return kFALSE;
+}
+
+Bool_t AliJetUnitArray::operator<( AliJetUnitArray unit) const
+{
+  // Less than operator used by sort
+  if( fUnitEnergy < unit.GetUnitEnergy())
+    return kTRUE;
+  else
+    return kFALSE;
+}
+
+Bool_t AliJetUnitArray::operator==( AliJetUnitArray unit) const
+{
+  // equality operator used by sort
+  if( fUnitEnergy == unit.GetUnitEnergy())
+    return kTRUE;
+  else
+    return kFALSE;
+}
diff --git a/JETAN/AliJetUnitArray.h b/JETAN/AliJetUnitArray.h
new file mode 100644 (file)
index 0000000..cc0f0d0
--- /dev/null
@@ -0,0 +1,101 @@
+#ifndef ALIUA1JETFINDERUNIT_H
+#define ALIUA1JETFINDERUNIT_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)
+
+#include <TObject.h>
+#include "AliJetFinderTypes.h"
+
+class AliJetUnitArray : public TObject
+{
+ public:
+  AliJetUnitArray();
+  ~AliJetUnitArray();
+
+  // Setter
+  void SetUnitEnergy(Float_t energy)  {fUnitEnergy = energy;}
+  void SetUnitEta(Float_t eta)        {fUnitEta = eta;} 
+  void SetUnitPhi(Float_t phi)        {fUnitPhi = phi;}
+  void SetUnitDeta(Float_t deta)      {fUnitDeta = deta;} 
+  void SetUnitDphi(Float_t dphi)      {fUnitDphi = dphi;}
+  void SetUnitID(Int_t id)            {fUnitID = id;}
+  void SetUnitEntries(Int_t num)      {fUnitNum = num;}
+  void SetUnitClusterID(Int_t id)     {fUnitClusterID = id;}
+  void SetUnitFlag(AliJetFinderUnitFlagType_t flag)    
+  { 
+    fUnitFlag = flag; 
+  } 
+  void SetUnitCutFlag(AliJetFinderUnitCutFlagType_t cutFlag)    
+  { 
+    fUnitCutFlag = cutFlag; 
+  } 
+  void SetUnitSignalFlag(AliJetFinderUnitSignalFlagType_t signalFlag)    
+  { 
+    fUnitSignalFlag = signalFlag; 
+  } 
+  void SetUnitDetectorFlag(AliJetFinderUnitDetectorFlagType_t detectorflag)    
+  { 
+    fUnitDetectorFlag = detectorflag; 
+  } 
+
+  // Getter
+  Float_t GetUnitEnergy() const            {return fUnitEnergy;}
+  Float_t GetUnitEta() const               {return fUnitEta;}
+  Float_t GetUnitPhi() const               {return fUnitPhi;}         
+  Float_t GetUnitDeta() const              {return fUnitDeta;}
+  Float_t GetUnitDphi() const              {return fUnitDphi;}         
+  Int_t   GetUnitID() const                {return fUnitID;}
+  Int_t   GetUnitEntries() const           {return fUnitNum;}
+  Int_t   GetUnitClusterID() const         {return fUnitClusterID;}
+  AliJetFinderUnitFlagType_t GetUnitFlag() const     
+  {
+         return fUnitFlag;
+  }
+  AliJetFinderUnitCutFlagType_t GetUnitCutFlag() const     
+  {
+         return fUnitCutFlag;
+  }
+  AliJetFinderUnitSignalFlagType_t GetUnitSignalFlag() const     
+  {
+         return fUnitSignalFlag;
+  }
+  AliJetFinderUnitDetectorFlagType_t GetUnitDetectorFlag() const     
+  {
+         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
+
+  ClassDef(AliJetUnitArray,3)
+};
+
+#endif
+