]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - JETAN/AliUA1JetFinder.cxx
- The part of JETAN dealing with ESD data has been separated from the one using MC...
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinder.cxx
index 2e291deaba34f97f986269cc32abec4750203ded..c19bb5d79d001cf64ded7b2d2bf189e608c69c7a 100755 (executable)
@@ -12,8 +12,6 @@
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-   
-  
  
 //---------------------------------------------------------------------
 // UA1 Jet finder 
@@ -24,6 +22,7 @@
 
 #include <Riostream.h>
 #include <TLorentzVector.h>
+#include <TFile.h>
 #include <TH2F.h>
 #include "AliUA1JetFinder.h"
 #include "AliUA1JetHeader.h"
@@ -37,14 +36,11 @@ ClassImp(AliUA1JetFinder)
 
 ////////////////////////////////////////////////////////////////////////
 
-AliUA1JetFinder::AliUA1JetFinder()
-
+AliUA1JetFinder::AliUA1JetFinder():
+  fHeader(0x0),
+  fLego(0x0)
 {
-  //
-  // Constructor
-  //
-  fHeader = 0x0;
-  fLego   = 0x0;
+  // Default constructor
 }
 
 ////////////////////////////////////////////////////////////////////////
@@ -52,13 +48,7 @@ AliUA1JetFinder::AliUA1JetFinder()
 AliUA1JetFinder::~AliUA1JetFinder()
 
 {
-  //
   // destructor
-  //
-
-  // delete fLego;            
-  
-  // reset and delete header
 }
 
 ////////////////////////////////////////////////////////////////////////
@@ -90,7 +80,6 @@ void AliUA1JetFinder::FindJets()
 {
   // initialize event, look for jets, download jet info
 
-
   // initialization in 2 steps
   // 1) transform input to pt,eta,phi plus lego
   // 2) dump lego
@@ -101,16 +90,19 @@ void AliUA1JetFinder::FindJets()
   if (nIn == 0) return;
     
   // local arrays for input
+  Float_t* enT  = new Float_t[nIn];
   Float_t* ptT  = new Float_t[nIn];
   Float_t* etaT = new Float_t[nIn];
   Float_t* phiT = new Float_t[nIn];
 
   // load input vectors
   for (Int_t i = 0; i < nIn; i++){
-      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
-      ptT[i]  = lv->Pt();
-      etaT[i] = lv->Eta();
-      phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
+    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());
+    if (fReader->GetCutFlag(i) == 1) 
       fLego->Fill(etaT[i], phiT[i], ptT[i]);
   }
   fJets->SetNinput(nIn);
@@ -125,8 +117,8 @@ void AliUA1JetFinder::FindJets()
   TAxis* xaxis = fLego->GetXaxis();
   TAxis* yaxis = fLego->GetYaxis();
   Float_t e = 0.0;
-  for (Int_t i = 1; i <= fHeader->GetNbinEta(); i++) {
-      for (Int_t j = 1; j <= fHeader->GetNbinPhi(); j++) {
+  for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
+      for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
          e = fLego->GetBinContent(i,j);
          if (e > 0.0) e -= fHeader->GetMinCellEt();
          if (e < 0.0) e = 0.;
@@ -140,7 +132,7 @@ void AliUA1JetFinder::FindJets()
   }
 
   // run the algo. Parameters from header
-  Int_t nTot      = (fHeader->GetNbinEta())*(fHeader->GetNbinPhi());
+  Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
   Float_t minmove = fHeader->GetMinMove();
   Float_t maxmove = fHeader->GetMaxMove();
   Int_t mode      = fHeader->GetMode();
@@ -149,60 +141,83 @@ void AliUA1JetFinder::FindJets()
   ua1_jet_finder(nCell, 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);
 
   // download jet info.   
-  AliJetReaderHeader* rheader = fReader->GetReaderHeader();
   for(Int_t i = 0; i < UA1JETS.njet; i++) {
     // reject events outside acceptable eta range
-    if (((UA1JETS.etaj[1][idx[i]])> (rheader->GetFiducialEtaMax()))
-       || ((UA1JETS.etaj[1][idx[i]]) < (rheader->GetFiducialEtaMin())))
-      continue;
-    Float_t px, py,pz,en; // 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);
-    fJets->AddJet(px, py, pz, en);
+    if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
+       || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
+         continue;
+      Float_t px, py,pz,en; // 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);
+      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[nIn];
-  for (Int_t i = 0; i < nIn; i++) injet[i]= -1;
+  Int_t* sflag = new Int_t[nIn];
+  for (Int_t i = 0; i < nIn; 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()];
 
   for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
-    Float_t pt_sig = 0.0;
-    mult[i] = 0;
-    for (Int_t j = 0; j < nIn; 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] == -1) {
-       if (fReader->GetSignalFlag(j) == 1) pt_sig+=ptT[j];
-       injet[j] = i;
-       mult[i]++;
+      Float_t pt_sig = 0.0;
+      mult[i] = 0;
+      ncell[i] = UA1JETS.ncellj[i];
+      for (Int_t j = 0; j < nIn; 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 ((fReader->GetCutFlag(j)) == 1 &&
+                 (etaT[j] < fHeader->GetLegoEtaMax()) &&
+                 (etaT[j] > fHeader->GetLegoEtaMin())) {
+                 injet[j] = i+1;
+                 mult[i]++;
+                 if (fReader->GetSignalFlag(j) == 1) {
+                   pt_sig+=ptT[j];
+                   sflag[j]=1;
+                 }
+             }
+         }
       }
-    }
-    percentage[i] = pt_sig/((Double_t) fJets->GetPt(i));    
+      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);
+  fJets->SetPhiIn(phiT);
+  fJets->SetPtIn(ptT);
+  fJets->SetEtAvg(UA1JETS.etavg);
+  delete ncell;
+  delete enT;
+  delete ptT;
+  delete etaT;
+  delete phiT;
+  delete injet;
+  delete idx;
+  delete mult;
+  delete percentage;
 }
 
 ////////////////////////////////////////////////////////////////////////
 
 void AliUA1JetFinder::Reset()
-
 {
   fLego->Reset();
   fJets->ClearJets();
@@ -216,17 +231,16 @@ void AliUA1JetFinder::WriteJHeaderToFile()
   fHeader->Write();
 }
 
-
 ////////////////////////////////////////////////////////////////////////
 
 void AliUA1JetFinder::Init()
 {
   // initializes some variables
   Float_t dEta, dPhi;
-  dEta = (fHeader->GetEtaMax()-fHeader->GetEtaMin())
-    /((Float_t) fHeader->GetNbinEta());
-  dPhi = (fHeader->GetPhiMax()-fHeader->GetPhiMin())
-    /((Float_t) fHeader->GetNbinPhi());
+  dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
+    /((Float_t) fHeader->GetLegoNbinEta());
+  dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
+    /((Float_t) fHeader->GetLegoNbinPhi());
 
   UA1CELL.etaCellSize = dEta;
   UA1CELL.phiCellSize = dPhi;
@@ -240,6 +254,7 @@ void AliUA1JetFinder::Init()
   // book lego
   fLego = new 
     TH2F("legoH","eta-phi",
-        fHeader->GetNbinEta(), fHeader->GetEtaMin(), fHeader->GetEtaMax(),
-        fHeader->GetNbinPhi(), fHeader->GetPhiMin(), fHeader->GetPhiMax());
+        fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(), 
+        fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(), 
+        fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
 }