New class for dE/dx analysis (comparison with Bethe-Bloch, check of response function...

[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinderV1.cxx

diff --git a/JETAN/AliUA1JetFinderV1.cxx b/JETAN/AliUA1JetFinderV1.cxx
index 4a625f0dd8ebfdf4a9b2e28f0dc8df2c19d0daea..6dafbdd15e411be2ab191561122ae69b547363ff 100644 (file)
--- a/JETAN/AliUA1JetFinderV1.cxx
+++ b/JETAN/AliUA1JetFinderV1.cxx
@@ -35,6 +35,7 @@
 #include "AliJetReader.h"
 #include "AliJet.h"
 #include "AliAODJet.h"
+#include "AliLog.h"
 
 
 ClassImp(AliUA1JetFinderV1)
@@ -161,11 +162,32 @@ void AliUA1JetFinderV1::FindJets()
      etbgTotal = etbgTotalN; // update with new background estimation
   } //end while
 
+  // add tracks to the jet if it wasn't yet done                                                         
+  if (header->GetBackgMode() == 0){
+    Float_t rc= header->GetRadius();
+    for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array                     
+      for(Int_t ijet=0; ijet<nj; ijet++){
+        Float_t deta = etaT[jpart] - etaJet[ijet];
+        Float_t dphi = phiT[jpart] - phiJet[ijet];
+       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+       if (dphi >  TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+        Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+       if(dr <= rc){ // particles inside this cone                                                    
+          injet[jpart] = ijet;
+          break;
+        }
+      }// end jets loop                                                                                
+    } //end particle loop                                                                              
+  }
+  
   // add jets to list
   Int_t* idxjets = new Int_t[nj];
   Int_t nselectj = 0;
 //  printf("Found %d jets \n", nj);
   
+  TRefArray *refs = 0;
+  Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
+  if (fromAod) refs = fReader->GetReferences();
   for(Int_t kj=0; kj<nj; kj++){
      if ((etaJet[kj] > (header->GetJetEtaMax())) ||
           (etaJet[kj] < (header->GetJetEtaMin())) ||
@@ -177,13 +199,21 @@ void AliUA1JetFinderV1::FindJets()
       en = TMath::Sqrt(px * px + py * py + pz * pz);
       fJets->AddJet(px, py, pz, en);
       AliAODJet jet(px, py, pz, en);
+
+      if (fromAod){
+        for(Int_t jpart = 0; jpart < nIn; jpart++) // loop for all particles in array
+          if (injet[jpart] == kj && fReader->GetCutFlag(jpart) == 1)
+                   jet.AddTrack(refs->At(jpart));  // check if the particle belongs to the jet and add the ref
+      }
+      
       //jet.Print("");
       
       AddJet(jet);
       
       idxjets[nselectj] = kj;
       nselectj++;
-  }
+  } //end particle loop
+
   //add signal percentage and total signal  in AliJets for analysis tool
   Float_t* percentage  = new Float_t[nselectj];
   Int_t* ncells      = new Int_t[nselectj];
@@ -245,21 +275,26 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t&
 {
 
    //dump lego
-  // check enough space! *to be done*
   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
-  const Int_t nBinsMax = 70000;
+  const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory
   
+  const Int_t nBinEta = header->GetLegoNbinEta();
+  const Int_t nBinPhi = header->GetLegoNbinPhi();
+  if((nBinPhi*nBinEta)>nBinsMax){
+    AliError("Too many bins of the ETA-PHI histogram");
+  }
+
   Float_t etCell[nBinsMax];   //! Cell Energy
   Float_t etaCell[nBinsMax];  //! Cell eta
   Float_t phiCell[nBinsMax];  //! Cell phi
-  Int_t   flagCell[nBinsMax]; //! Cell flag
+  Short_t   flagCell[nBinsMax]; //! Cell flag
 
   Int_t nCell = 0;
   TAxis* xaxis = fLego->GetXaxis();
   TAxis* yaxis = fLego->GetYaxis();
   Float_t e = 0.0;
-  for (Int_t i = 1; i <= header->GetLegoNbinEta(); i++) {
-      for (Int_t j = 1; j <= header->GetLegoNbinPhi(); j++) {
+  for (Int_t i = 1; i <= nBinEta; i++) {
+      for (Int_t j = 1; j <= nBinPhi; j++) {
               e = fLego->GetBinContent(i,j);
               if (e < 0.0) continue; // don't include this cells
               Float_t eta  = xaxis->GetBinCenter(i);