Build list of tracks associated to jet.

diff --git a/EMCAL/AliEMCALJetFinder.cxx b/EMCAL/AliEMCALJetFinder.cxx
index 45da980d08a4f8b3a5899f7adac5de765e03cfd4..b5ff17a2cd9b423c047c6d06a71b20c25e8b6e10 100644 (file)
--- a/EMCAL/AliEMCALJetFinder.cxx
+++ b/EMCAL/AliEMCALJetFinder.cxx
@@ -36,9 +36,13 @@ ClassImp(AliEMCALJetFinder)
 AliEMCALJetFinder::AliEMCALJetFinder()
 {
 // Default constructor
-    fJets  = 0;
-    fNjets = 0;
-    fLego  = 0;
+    fJets      = 0;
+    fNjets     = 0;
+    fLego      = 0;
+    fTrackList = 0;
+    fPtT       = 0;
+    fEtaT      = 0;
+    fPhiT      = 0;
 }
 
 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
@@ -54,6 +58,11 @@ AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
         fPhiCell[i] = 0.;
     }
     fLego = 0;
+    fTrackList = 0;
+    fTrackList = 0;
+    fPtT       = 0;
+    fEtaT      = 0;
+    fPhiT      = 0;
 }
 
 
@@ -72,6 +81,7 @@ AliEMCALJetFinder::~AliEMCALJetFinder()
 
 
 
+
 #ifndef WIN32
 # define jet_finder_ua1 jet_finder_ua1_
 # define hf1 hf1_
@@ -123,6 +133,14 @@ void AliEMCALJetFinder::Find()
     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell, 
                   min_move, max_move, mode, prec_bg, ierror);
     // Write result to output
+    Int_t njet = Njets();
+    for (Int_t nj=0; nj<njet; nj++)
+    {
+       fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
+                                   JetPhiW(nj),
+                                   JetEtaW(nj));
+    }
+    FindTracksInJetCone();
     WriteJets();
 }
 
@@ -175,6 +193,7 @@ void AliEMCALJetFinder::SetConeRadius(Float_t par)
 {
 // Set jet cone radius
     EMCALJETPARAM.coneRad = par;
+    fConeRadius = par;
 }
 
 void AliEMCALJetFinder::SetEtSeed(Float_t par)
@@ -195,6 +214,12 @@ void AliEMCALJetFinder::SetMinCellEt(Float_t par)
     EMCALJETPARAM.etMin = par;
 }
 
+void AliEMCALJetFinder::SetPtCut(Float_t par)
+{
+// Set pT cut on charged tracks
+    fPtCut = par;
+}
+
 
 void AliEMCALJetFinder::Test()
 {
@@ -261,14 +286,12 @@ void AliEMCALJetFinder::WriteJets()
                                 file);
     }
     Int_t njet = Njets();
-    for (Int_t nj=0; nj<njet; nj++)
+    for (Int_t nj = 0; nj < njet; nj++)
     {
-       AliEMCALJet* jet = new AliEMCALJet(JetEnergy(nj),
-                                          JetPhiW(nj),
-                                          JetEtaW(nj));
-       AddJet(*jet);
-       delete jet;
+       AddJet(*fJetT[nj]);
+       delete fJetT[nj];
     }
+
     Int_t nev = gAlice->GetHeader()->GetEvent();
     gAlice->TreeR()->Fill();
     char hname[30];
@@ -352,21 +375,43 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 //
 // Access particle information    
     Int_t npart = (gAlice->GetHeader())->GetNprimary();
-    for (Int_t part=2; part<npart; part++) {
+// Create track list
+//
+// 0: not selected
+// 1: selected for jet finding
+// 2: inside jet
+// ....
+    if (fTrackList) delete[] fTrackList;
+    if (fPtT)       delete[] fPtT;
+    if (fEtaT)      delete[] fEtaT;
+    if (fPhiT)      delete[] fPhiT;
+    
+    fTrackList = new Int_t  [npart];
+    fPtT       = new Float_t[npart];
+    fEtaT      = new Float_t[npart];
+    fPhiT      = new Float_t[npart];
+    fNt        = npart;
+
+    for (Int_t part = 0; part < npart; part++) {
        TParticle *MPart = gAlice->Particle(part);
        Int_t mpart   = MPart->GetPdgCode();
        Int_t child1  = MPart->GetFirstDaughter();
        Float_t pT    = MPart->Pt();
        Float_t phi   = MPart->Phi();
-       Float_t theta = MPart->Theta();
-       Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
+       Float_t eta   = MPart->Eta();
 //     if (part == 6 || part == 7)
 //     {
 //         printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f", 
 //                part-5, pT, eta, phi);
 //     }
-           
-       if (pT == 0.) continue;
+       
+       fTrackList[part] = 0;
+       fPtT[part]       = pT;
+       fEtaT[part]      = eta;
+       fPhiT[part]      = phi;
+
+       if (part < 2) continue;
+       if (pT == 0 || pT < fPtCut) continue;
 // charged or neutral 
        if (ich == 0) {
            TParticlePDG* pdgP = MPart->GetPDG();
@@ -383,6 +428,7 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
        if (child1 >= 0 && child1 < npart) continue;
 //     printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f", 
 //     part, mpart, child1, eta, phi, pT);
+       fTrackList[part] = 1;
        fLego->Fill(eta, phi, pT);
     } // primary loop
     DumpLego();
@@ -430,10 +476,10 @@ void AliEMCALJetFinder::FillFromHits(Int_t flag)
 //         Int_t   index  =    mHit->GetId();     // cell index
 //         Float_t eta, phi;
 //         geom->EtaPhiFromIndex(index,  eta, phi);
-           Float_t r     = TMath::Sqrt(x*x+y*y);
-           Float_t theta = TMath::ATan2(r,z);
-           Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
-           Float_t phi   = TMath::ATan2(y,x);
+           Float_t r      = TMath::Sqrt(x*x+y*y);
+           Float_t theta  = TMath::ATan2(r,z);
+           Float_t eta    = -TMath::Log(TMath::Tan(theta/2.));
+           Float_t phi    = TMath::ATan2(y,x);
            fLego->Fill(eta, phi, eloss);
 //         if (eloss > 1.) printf("\nx,y,z %f %f %f %f %f", 
 //         r, z, eta, phi, eloss);
@@ -443,6 +489,51 @@ void AliEMCALJetFinder::FillFromHits(Int_t flag)
     DumpLego();
 }
 
+void AliEMCALJetFinder::FindTracksInJetCone()
+{
+//
+//  Build list of tracks inside jet cone
+//
+//  Loop over jets
+    Int_t njet = Njets();
+    for (Int_t nj = 0; nj < njet; nj++)
+    {
+       Float_t etaj = JetEtaW(nj);
+       Float_t phij = JetPhiW(nj);     
+       Int_t   nT   = 0;           // number of associated tracks
+       
+// Loop over particles
+       for (Int_t part = 0; part < fNt; part++) {
+           if (!fTrackList[part]) continue;
+           Float_t phi      = fPhiT[part];
+           Float_t eta      = fEtaT[part];
+           Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
+                                          (phij-phi)*(phij-phi));
+           if (dr < fConeRadius) {
+               fTrackList[part] = nj+2;
+               nT++;
+           } // < ConeRadius ?
+       } // particle loop
+       Float_t* ptT  = new Float_t[nT];
+       Float_t* etaT = new Float_t[nT];
+       Float_t* phiT = new Float_t[nT];
+       Int_t    iT   = 0;
+       for (Int_t part = 0; part < fNt; part++) {
+           if (fTrackList[part] == nj+2) {
+               ptT [iT] = fPtT [part];
+               etaT[iT] = fEtaT[part];
+               phiT[iT] = fPhiT[part];
+               iT++;
+           } // particle associated
+       } // particle loop
+       fJetT[nj]->SetTrackList(nT, ptT, etaT, phiT);
+       delete[] ptT;
+       delete[] etaT;
+       delete[] phiT;
+    } // jet loop loop
+}
+
+
 
 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
 {

diff --git a/EMCAL/AliEMCALJetFinder.h b/EMCAL/AliEMCALJetFinder.h
index e8d95267399c02cf0c7876aa99c84027000bc212..7ee67af6e4030cae4f616b4fb81d5df0085806eb 100644 (file)
--- a/EMCAL/AliEMCALJetFinder.h
+++ b/EMCAL/AliEMCALJetFinder.h
@@ -16,56 +16,67 @@ class TH2F;
 
 class AliEMCALJetFinder : public TTask {
  public:
-  AliEMCALJetFinder();
-  AliEMCALJetFinder(const char* name, const char *title);
-  virtual ~AliEMCALJetFinder();
-  virtual void Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000], 
-                   Float_t etac[30000], Float_t phic[30000],
-                   Float_t min_move, Float_t max_move, Int_t mode,
-                   Float_t prec_bg, Int_t ierror);
-  virtual void Find();
-  virtual void Test();
-  // Geometry
-  virtual void SetCellSize(Float_t eta, Float_t phi);
-  // Parameters
-  virtual void SetConeRadius(Float_t par);
-  virtual void SetEtSeed(Float_t par);
-  virtual void SetMinJetEt(Float_t par);
-  virtual void SetMinCellEt(Float_t par);
-  // Results
-  virtual Int_t   Njets();
-  virtual Float_t JetEnergy(Int_t);
-  virtual Float_t JetPhiL(Int_t);
-  virtual Float_t JetPhiW(Int_t);
-  virtual Float_t JetEtaL(Int_t);  
-  virtual Float_t JetEtaW(Int_t);
-  virtual TH2F* GetLego() {return fLego;}
-  // I/O
-  virtual void FillFromHits(Int_t flag = 0);
-  virtual void FillFromTracks(Int_t flag = 0, Int_t ich = 0);
-  virtual void AddJet(const AliEMCALJet& jet);
-  virtual void WriteJets();
-  virtual void ResetJets();
-  virtual TClonesArray*  Jets() {return fJets;}
+    AliEMCALJetFinder();
+    AliEMCALJetFinder(const char* name, const char *title);
+    virtual ~AliEMCALJetFinder();
+    virtual void Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000], 
+                     Float_t etac[30000], Float_t phic[30000],
+                     Float_t min_move, Float_t max_move, Int_t mode,
+                     Float_t prec_bg, Int_t ierror);
+    virtual void Find();
+    virtual void FindTracksInJetCone();
+    virtual void Test();
+    // Geometry
+    virtual void SetCellSize(Float_t eta, Float_t phi);
+    // Parameters
+    virtual void SetConeRadius(Float_t par);
+    virtual void SetEtSeed(Float_t par);
+    virtual void SetMinJetEt(Float_t par);
+    virtual void SetMinCellEt(Float_t par);
+    virtual void SetPtCut(Float_t par);
+    // Results
+    virtual Int_t   Njets();
+    virtual Float_t JetEnergy(Int_t);
+    virtual Float_t JetPhiL(Int_t);
+    virtual Float_t JetPhiW(Int_t);
+    virtual Float_t JetEtaL(Int_t);  
+    virtual Float_t JetEtaW(Int_t);
+    virtual TH2F* GetLego() {return fLego;}
+    // I/O
+    virtual void FillFromHits(Int_t flag = 0);
+    virtual void FillFromTracks(Int_t flag = 0, Int_t ich = 0);
+    virtual void AddJet(const AliEMCALJet& jet);
+    virtual void WriteJets();
+    virtual void ResetJets();
+    virtual TClonesArray* Jets() {return fJets;}
  private:
-  virtual void BookLego();
-  virtual void DumpLego();
-  virtual void ResetMap();
+    virtual void BookLego();
+    virtual void DumpLego();
+    virtual void ResetMap();
  protected:
-  TClonesArray*      fJets;            //! List of Jets
-  TH2F*              fLego;            //! Lego Histo 
-  Int_t              fNjets;           //! Number of Jets
-  Float_t            fDeta;            //! eta cell size 
-  Float_t            fDphi;            //! phi cell size
-  Int_t              fNcell;           //! number of cells
-  Int_t              fNtot;            //! total number of cells
-  Int_t              fNbinEta;         //! number of cells in eta
-  Int_t              fNbinPhi;         //! number of cells in phi
-  Float_t            fEtCell[30000];   //! Cell Energy
-  Float_t            fEtaCell[30000];  //! Cell eta
-  Float_t            fPhiCell[30000];  //! Cell phi
-  
-  ClassDef(AliEMCALJetFinder,1) // JetFinder for EMCAL
-} ;
+    TClonesArray*      fJets;            //! List of Jets
+    TH2F*              fLego;            //! Lego Histo
+    AliEMCALJet*       fJetT[10];        //! Jet temporary storage
+    Float_t            fConeRadius;      //  Cone radius
+    Float_t            fPtCut;           //  Pt cut on charged tracks 
+    Int_t              fNjets;           //! Number of Jets
+    Float_t            fDeta;            //! eta cell size 
+    Float_t            fDphi;            //! phi cell size
+    Int_t              fNcell;           //! number of cells
+    Int_t              fNtot;            //! total number of cells
+    Int_t              fNbinEta;         //! number of cells in eta
+    Int_t              fNbinPhi;         //! number of cells in phi
+    Float_t            fEtCell[30000];   //! Cell Energy
+    Float_t            fEtaCell[30000];  //! Cell eta
+    Float_t            fPhiCell[30000];  //! Cell phi
+    Int_t*             fTrackList;       //! List of selected tracks
+    Int_t              fNt;              //! number of tracks
+    Float_t*           fPtT;             //! Pt   of tracks in jet cone
+    Float_t*           fEtaT;            //! Eta  of tracks in jet cone
+    Float_t*           fPhiT;            //! Phi  of tracks in jet cone
+    ClassDef(AliEMCALJetFinder,2)        // JetFinder for EMCAL
+}
+;
 #endif // ALIEMCALJetFinder_H
 
+