]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - JETAN/AliFastJetFinder.cxx
AddTask fix
[u/mrichter/AliRoot.git] / JETAN / AliFastJetFinder.cxx
index 833bfa5815e4e7d1125bfa2c144a54f9a8bdd5b4..f15989bffd59ea8ea7adb715a2276608bcf8ac6b 100644 (file)
 
 //---------------------------------------------------------------------
 // FastJet v2.3.4 finder algorithm interface
+// Last modification: Neutral cell energy included in the jet reconstruction
+//
+// Authors: Rafael.Diaz.Valdes@cern.ch
+//          Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option)
 //
-// Author: Rafael.Diaz.Valdes@cern.ch
-//  
 //---------------------------------------------------------------------
 
+
 #include <Riostream.h>
 #include <TLorentzVector.h>
 #include <TFile.h>
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TArrayF.h>
-#include <TRandom.h>
 #include <TClonesArray.h>
 
 #include "AliFastJetFinder.h"
-#include "AliFastJetHeader.h"
+#include "AliFastJetHeaderV1.h"
 #include "AliJetReaderHeader.h"
 #include "AliJetReader.h"
-#include "AliJet.h"
-
+#include "AliJetUnitArray.h"
+#include "AliFastJetInput.h"
+#include "AliJetBkg.h"
+#include "AliAODJetEventBackground.h"
+#include "AliAODTrack.h"
 
 #include "fastjet/PseudoJet.hh"
 #include "fastjet/ClusterSequenceArea.hh"
@@ -56,10 +61,14 @@ using namespace std;
 
 
 ClassImp(AliFastJetFinder)
+
+
 //____________________________________________________________________________
 
 AliFastJetFinder::AliFastJetFinder():
-    AliJetFinder()
+  AliJetFinder(),
+  fInputFJ(new AliFastJetInput()),
+  fJetBkg(new  AliJetBkg())
 {
   // Constructor
 }
@@ -69,132 +78,222 @@ AliFastJetFinder::AliFastJetFinder():
 AliFastJetFinder::~AliFastJetFinder()
 {
   // destructor
+  delete  fInputFJ; fInputFJ = 0;
+  delete  fJetBkg; fJetBkg = 0;
 }
 
 //______________________________________________________________________________
 void AliFastJetFinder::FindJets()
 {
-  
-  Bool_t debug = kFALSE;
-  
+
   //pick up fastjet header
-  AliFastJetHeader *header = (AliFastJetHeader*)fHeader;
+  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+  Int_t  debug  = header->GetDebug();     // debug option
+  Bool_t bgMode = header->GetBGMode();    // choose to subtract BG or not
+  if(debug)cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl;
 
   // check if we are reading AOD jets
   TRefArray *refs = 0;
   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
   if (fromAod) { refs = fReader->GetReferences(); }
-  
+
   // RUN ALGORITHM  
   // read input particles -----------------------------
-  vector<fastjet::PseudoJet> input_particles;
-  TClonesArray *lvArray = fReader->GetMomentumArray();
-  if(lvArray == 0) { cout << "Could not get the momentum array" << endl; return; }
-  Int_t nIn =  lvArray->GetEntries();
-  if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }
-  //Int_t nJets = 0; // n jets in this event
-  fJets->SetNinput(nIn) ; // number of input objects
-  Float_t px,py,pz,en;
-  // load input vectors
-  for(Int_t i = 0; i < nIn; i++){ // loop for all input particles
-      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
-      px = lv->Px();
-      py = lv->Py();
-      pz = lv->Pz();
-      en = lv->Energy();
-      
-      fastjet::PseudoJet input_part(px,py,pz,en); // create PseudoJet object
-      input_part.set_user_index(i); //label the particle into Fastjet algortihm
-      input_particles.push_back(input_part);  // back of the input_particles vector  
-  } // end loop 
-  
-  
+
+  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
+  if(inputParticles.size()==0){
+    if(debug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
+    return;
+  }
+
   // create an object that represents your choice of jet algorithm, and 
   // the associated parameters
-  double Rparam = header->GetRparam();
+  double rParam = header->GetRparam();
+  double rBkgParam = header->GetRparamBkg();
   fastjet::Strategy strategy = header->GetStrategy();
-  fastjet::RecombinationScheme recomb_scheme = header->GetRecombScheme();
+  fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
   fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); 
-  fastjet::JetDefinition jet_def(algorithm, Rparam, recomb_scheme, strategy);
-  
+  fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy);
+
   // create an object that specifies how we to define the area
-  fastjet::AreaDefinition area_def;
-  double ghost_etamax = header->GetGhostEtaMax(); 
-  double ghost_area   = header->GetGhostArea(); 
-  int    active_area_repeats = header->GetActiveAreaRepeats(); 
+  fastjet::AreaDefinition areaDef;
+  double ghostEtamax = header->GetGhostEtaMax(); 
+  double ghostArea   = header->GetGhostArea(); 
+  int    activeAreaRepeats = header->GetActiveAreaRepeats(); 
   
   // now create the object that holds info about ghosts
-  fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, ghost_area);
+  fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
   // and from that get an area definition
-  fastjet::AreaType area_type = header->GetAreaType();
-  area_def = fastjet::AreaDefinition(area_type,ghost_spec);
+  fastjet::AreaType areaType = header->GetAreaType();
+  areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
   
+  //***************************** JETS FINDING
   // run the jet clustering with the above jet definition
-  fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def);
-  
-  
-  // save a comment in the header
-  TString comment = "Running FastJet algorithm with the following setup. ";
-  comment+= "Jet definition: ";
-  comment+= TString(jet_def.description());
-  comment+= ". Area definition: ";
-  comment+= TString(area_def.description());
-  comment+= ". Strategy adopted by FastJet: ";
-  comment+= TString(clust_seq.strategy_string());
-  header->SetComment(comment);
-  if(debug){
-    cout << "--------------------------------------------------------" << endl;
-    cout << comment << endl;
-    cout << "--------------------------------------------------------" << endl;
+  fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
+
+  vector<fastjet::PseudoJet> jets;
+
+  if(bgMode) // Do BG subtraction directly with the same algorithm (cambridge or kt) for jet signal and background
+    {
+      //***************************** CLUSTER JETS FINDING FOR RHO ESTIMATION
+      // run the jet clustering with the above jet definition
+         fastjet::JetAlgorithm algorithmBkg = header->GetBGAlgorithm();
+         fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy);
+      fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef);
+
+      // save a comment in the header
+      TString comment = "Running FastJet algorithm with the following setup. ";
+      comment+= "Jet definition: ";
+      comment+= TString(jetDef.description());
+      comment+= "Jet bckg definition: ";
+      comment+= TString(jetDefBkg.description());
+      comment+= ". Area definition: ";
+      comment+= TString(areaDef.description());
+      comment+= ". Strategy adopted by FastJet and bkg: ";
+      comment+= TString(clust_seq.strategy_string());
+      header->SetComment(comment);
+      if(debug){
+          cout << "--------------------------------------------------------" << endl;
+          cout << comment << endl;
+          cout << "--------------------------------------------------------" << endl;
+      }
+      //header->PrintParameters();
+      
+      // extract the inclusive jets with pt > ptmin, sorted by pt
+      double ptmin = header->GetPtMin(); 
+      vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets();
+
+      //subtract background // ===========================================
+      // set the rapididty , phi range within which to study the background 
+      double rapMax = header->GetRapMax(); 
+      double rapMin = header->GetRapMin();
+      double phiMax = header->GetPhiMax();
+      double phiMin = header->GetPhiMin();
+      fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
+      
+      // Extract rho and sigma
+      Double_t rho = 0.;
+      Double_t sigma = 0.;
+      Double_t meanarea = 0.;
+      Bool_t Use4VectorArea = header->Use4VectorArea();
+      vector<fastjet::PseudoJet> bkgJets = clust_seq_bkg.inclusive_jets();
+      clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, Use4VectorArea, rho, sigma, meanarea, false);
+
+         // subtract background and extract jets bkg subtracted
+      vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptmin);
+
+         // print out
+      //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n";
+      //cout << "---------------------------------------\n";
+      //cout << endl;
+      //printf(" ijet   rap      phi        Pt         area  +-   err\n");
+      
+      // sort jets into increasing pt
+      jets = sorted_by_pt(subJets);  
+
+    }
+  else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
+
+         // save a comment in the header
+         TString comment = "Running FastJet algorithm with the following setup. ";
+         comment+= "Jet definition: ";
+         comment+= TString(jetDef.description());
+         comment+= ". Strategy adopted by FastJet: ";
+         comment+= TString(clust_seq.strategy_string());
+         header->SetComment(comment);
+         if(debug){
+                 cout << "--------------------------------------------------------" << endl;
+                 cout << comment << endl;
+                 cout << "--------------------------------------------------------" << endl;
+         }
+         //header->PrintParameters();
+
+         // extract the inclusive jets with pt > ptmin, sorted by pt
+         double ptmin = header->GetPtMin();
+         vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
+
+         jets = sorted_by_pt(inclusiveJets); // Added by me
+
   }
-  //header->PrintParameters();
-  
-  
-  // extract the inclusive jets with pt > ptmin, sorted by pt
-  double ptmin = header->GetPtMin(); 
-  vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
-  
-  //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
-  //subtract background // ===========================================
-  // set the rapididty , phi range within which to study the background 
-  double rap_max = header->GetRapMax(); 
-  double rap_min = header->GetRapMin();
-  double phi_max = header->GetPhiMax();
-  double phi_min = header->GetPhiMin();
-  fastjet::RangeDefinition range(rap_min, rap_max, phi_min, phi_max);
-  // subtract background
-  vector<fastjet::PseudoJet> sub_jets =  clust_seq.subtracted_jets(range,ptmin);  
-  
-  // print out
-  //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n";
-  //cout << "---------------------------------------\n";
-  //cout << endl;
-  //printf(" ijet   rap      phi        Pt         area  +-   err\n");
-   
-  // sort jets into increasing pt
-  vector<fastjet::PseudoJet> jets = sorted_by_pt(sub_jets);  
+
   for (size_t j = 0; j < jets.size(); j++) { // loop for jets
 
-//    double area     = clust_seq.area(jets[j]);
-//    double area_error = clust_seq.area_error(jets[j]);
-
-    //printf("%5u %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, area_error);
-       
-       // go to write AOD  info
-    AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
-   //cout << "Printing jet " << endl;
-    if(debug) aodjet.Print("");
-   //cout << "Adding jet ... " ;
-    AddJet(aodjet);
-   //cout << "added \n" << endl;
-   
-  } // end loop for jets
+         double area      = clust_seq.area(jets[j]);
+         double areaError = clust_seq.area_error(jets[j]);
+
+         if(debug) printf("Jet found %5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, areaError);
+
+      vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
+      int nCon= constituents.size();
+      TArrayI ind(nCon);
+      
+      if ((jets[j].eta() > (header->GetJetEtaMax())) ||
+          (jets[j].eta() < (header->GetJetEtaMin())) ||
+          (jets[j].phi() > (header->GetJetPhiMax())) ||
+          (jets[j].phi() < (header->GetJetPhiMin())) ||
+          (jets[j].perp() < header->GetPtMin())) continue; // acceptance eta range and etmin
+
+         // go to write AOD  info
+         AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
+      aodjet.SetEffArea(area,areaError);
+         //cout << "Printing jet " << endl;
+         if(debug) aodjet.Print("");
+
+         Int_t count=0;
+         for (int i=0; i < nCon; i++)
+         {
+                 fastjet::PseudoJet mPart=constituents[i];
+                 ind[i]=mPart.user_index();
+                 //    cout<<i<<"  index="<<ind[i]<<endl;
+
+                 if(fromAod)
+                 {
+                         if(fReader->GetReaderHeader()->GetDetector()==0)
+                         {
+                                 for (Int_t iref = 0; iref < refs->GetEntries(); iref++)
+                                 {
+                                         if(iref==ind[i]){
+                                                 AliAODTrack * track = (AliAODTrack*)refs->At(iref);
+                                                 aodjet.AddTrack(track);
+                                         }
+                                 }
+
+                         } else {
+
+                                 TClonesArray* fUnit = fReader->GetUnitArray(); //Big mmentum array
+                                 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
+                                 Int_t         nIn = fUnit->GetEntries();
+
+                                 //internal loop over all the unit cells
+                                 Int_t ipart = 0;
+
+                                 for(Int_t ii=0; ii<nIn; ii++)
+                                 {
+                                         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(ii);
+                                         if(uArray->GetUnitEnergy()>0.){
+                                                 uArray->SetUnitTrackID(ipart);//used to have the index available in AliJEtBkg
+                                                 if(ipart==ind[i]){
+                                                         TRefArray* trackArray = (TRefArray*)uArray->GetUnitTrackRef();
+                                                         Int_t tracksInCell = trackArray->GetEntries();
+                                                         for(int ji = 0; ji < tracksInCell; ji++){
+                                                                 AliAODTrack * track = (AliAODTrack*)trackArray->At(ji);
+                                                                 aodjet.AddTrack(track);
+                                                         }
+
+                                                         count++;
+                                                 }
+                                                 ipart++;
+                                         }
+                                 }
+                         }
+                 }
+        } 
+
+ AddJet(aodjet);
+
+
+  } // End loop on jets
 
-    
 }
 
 //____________________________________________________________________________
@@ -204,7 +303,7 @@ void AliFastJetFinder::RunTest(const char* datafile)
 
    // This simple test run the kt algorithm for an ascii file testdata.dat
    // read input particles -----------------------------
-  vector<fastjet::PseudoJet> input_particles;
+  vector<fastjet::PseudoJet> inputParticles;
   Float_t px,py,pz,en;
   ifstream in;
   Int_t nlines = 0;
@@ -216,7 +315,7 @@ void AliFastJetFinder::RunTest(const char* datafile)
       if (!in.good()) break;
       //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz);
       nlines++;
-      input_particles.push_back(fastjet::PseudoJet(px,py,pz,en)); 
+      inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en)); 
    }
    //printf(" found %d pointsn",nlines);
    in.close();
@@ -224,66 +323,66 @@ void AliFastJetFinder::RunTest(const char* datafile)
  
   // create an object that represents your choice of jet algorithm, and 
   // the associated parameters
-  double Rparam = 1.0;
+  double rParam = 1.0;
   fastjet::Strategy strategy = fastjet::Best;
-  fastjet::RecombinationScheme recomb_scheme = fastjet::BIpt_scheme;
-  fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, recomb_scheme, strategy);
+  fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme;
+  fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy);
   
   
  
   // create an object that specifies how we to define the area
-  fastjet::AreaDefinition area_def;
-  double ghost_etamax = 7.0;
-  double ghost_area    = 0.05;
-  int    active_area_repeats = 1;
+  fastjet::AreaDefinition areaDef;
+  double ghostEtamax = 7.0;
+  double ghostArea    = 0.05;
+  int    activeAreaRepeats = 1;
   
 
   // now create the object that holds info about ghosts
-  fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, ghost_area);
+  fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
   // and from that get an area definition
-  area_def = fastjet::AreaDefinition(fastjet::active_area,ghost_spec);
+  areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec);
   
 
   // run the jet clustering with the above jet definition
-  fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def);
+  fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
   
   
   // tell the user what was done
   cout << "--------------------------------------------------------" << endl;
-  cout << "Jet definition was: " << jet_def.description() << endl;
-  cout << "Area definition was: " << area_def.description() << endl;
+  cout << "Jet definition was: " << jetDef.description() << endl;
+  cout << "Area definition was: " << areaDef.description() << endl;
   cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl;
   cout << "--------------------------------------------------------" << endl;
  
   
   // extract the inclusive jets with pt > 5 GeV, sorted by pt
   double ptmin = 5.0;
-  vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
+  vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
   
   cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
  
  
   //subtract background // ===========================================
   // set the rapididty range within which to study the background 
-  double rap_max = ghost_etamax - Rparam;
-  fastjet::RangeDefinition range(rap_max);
+  double rapMax = ghostEtamax - rParam;
+  fastjet::RangeDefinition range(rapMax);
   // subtract background
-  vector<fastjet::PseudoJet> sub_jets =  clust_seq.subtracted_jets(range,ptmin);  
+  vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);  
   
   // print them out //================================================
   cout << "Printing inclusive jets  after background subtraction \n";
   cout << "------------------------------------------------------\n";
   // sort jets into increasing pt
-  vector<fastjet::PseudoJet> jets = sorted_by_pt(sub_jets);  
+  vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);  
 
   printf(" ijet   rap      phi        Pt         area  +-   err\n");
   for (size_t j = 0; j < jets.size(); j++) {
 
     double area     = clust_seq.area(jets[j]);
-    double area_error = clust_seq.area_error(jets[j]);
+    double areaError = clust_seq.area_error(jets[j]);
 
-    printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",j,jets[j].rap(),
-          jets[j].phi(),jets[j].perp(), area, area_error);
+    printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(),
+          jets[j].phi(),jets[j].perp(), area, areaError);
   }
   cout << endl;
   // ================================================================
@@ -294,9 +393,152 @@ void AliFastJetFinder::RunTest(const char* datafile)
 
 //____________________________________________________________________________
 
-void AliFastJetFinder::WriteJHeaderToFile()
+void AliFastJetFinder::WriteJHeaderToFile() const
 {
   fHeader->Write();
 }
 
 //____________________________________________________________________________
+
+Float_t  AliFastJetFinder::EtaToTheta(Float_t arg)
+{
+  //  return (180./TMath::Pi())*2.*atan(exp(-arg));
+  return 2.*atan(exp(-arg));
+
+
+}
+
+//____________________________________________________________________________
+
+void AliFastJetFinder::InitTask(TChain *tree)
+{
+
+  printf("Fast jet finder initialization ******************");
+  fReader->CreateTasks(tree);
+
+}
+
+
+Bool_t AliFastJetFinder::ProcessEvent()
+{
+  //
+  // Process one event
+  // from meomntum array
+
+  Bool_t ok = fReader->FillMomentumArray();
+
+  if (!ok) return kFALSE;
+  fInputFJ->SetHeader(fHeader);
+  fInputFJ->SetReader(fReader);
+  fInputFJ->FillInput();
+  // Jets
+  FindJets(); 
+
+
+  if( fAODEvBkg){
+    fJetBkg->SetHeader(fHeader);
+    fJetBkg->SetReader(fReader);
+    Double_t sigma1 = 0,meanarea1= 0,sigma2 = 0,meanarea2 = 0;
+    Double_t bkg1 = 0,bkg2 = 0;
+    
+    fJetBkg->SetFastJetInput(fInputFJ);
+    fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
+    fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
+    fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+    fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2);
+  }
+
+
+
+
+  /*
+  fJetBkg->SetFastJetInput(fInputFJ);
+  Double_t bkg1=fJetBkg->BkgFastJet();
+  Double_t bkg2=fJetBkg->BkgChargedFastJet();
+  Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets);
+  Double_t bkg4=fJetBkg->BkgRemoveJetLeading(fAODjets);
+  
+  fAODEvBkg->SetBackground(0,bkg1);
+  fAODEvBkg->SetBackground(1,bkg2);
+  fAODEvBkg->SetBackground(2,bkg3);
+  fAODEvBkg->SetBackground(3,bkg4);
+  */
+  Reset();  
+  return kTRUE;
+
+}
+
+Bool_t AliFastJetFinder::ProcessEvent2()
+{
+  //
+  // Process one event
+  // Charged only or charged+neutral jets
+  //
+
+  TRefArray* ref = new TRefArray();
+  Bool_t procid = kFALSE;
+  Bool_t ok = fReader->ExecTasks(procid,ref);
+
+  // Delete reference pointer  
+  if (!ok) {delete ref; return kFALSE;}
+  
+  // Leading particles
+  fInputFJ->SetHeader(fHeader);
+  fInputFJ->SetReader(fReader);
+  fInputFJ->FillInput();
+  
+  // Jets
+  FindJets();
+  
+  if( fAODEvBkg){
+    fJetBkg->SetHeader(fHeader);
+    fJetBkg->SetReader(fReader);
+    fJetBkg->SetFastJetInput(fInputFJ);
+    Double_t sigma1,meanarea1,sigma2,meanarea2;
+    Double_t bkg1,bkg2;
+    fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
+    fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
+    fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+    fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2);
+  }
+
+
+//  Double_t bkg1=fJetBkg->BkgFastJet();
+//  Double_t bkg2=fJetBkg->BkgChargedFastJet();
+//  Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets);
+//  Double_t bkg4=fJetBkg->BkgRemoveJetLeading(fAODjets);
+  
+//  fAODEvBkg->SetBackground(0,bkg1);
+//  fAODEvBkg->SetBackground(1,bkg2);
+//  fAODEvBkg->SetBackground(2,bkg3);
+//  fAODEvBkg->SetBackground(3,bkg4);
+  
+  Int_t nEntRef    = ref->GetEntries();
+
+  for(Int_t i=0; i<nEntRef; i++)
+    { 
+      // Reset the UnitArray content which were referenced
+      ((AliJetUnitArray*)ref->At(i))->SetUnitTrackID(0);
+      ((AliJetUnitArray*)ref->At(i))->SetUnitEnergy(0.);
+      ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag(kPtSmaller);
+      ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag2(kPtSmaller);
+      ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlag(kBad);
+      ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlagC(kTRUE,kBad);
+      ((AliJetUnitArray*)ref->At(i))->SetUnitDetectorFlag(kTpc);
+      ((AliJetUnitArray*)ref->At(i))->SetUnitFlag(kOutJet);
+      ((AliJetUnitArray*)ref->At(i))->ClearUnitTrackRef();
+
+      // Reset process ID
+      AliJetUnitArray* uA = (AliJetUnitArray*)ref->At(i);
+      uA->ResetBit(kIsReferenced);
+      uA->SetUniqueID(0);     
+    }
+
+  // Delete the reference pointer
+  ref->Delete();
+  delete ref;
+
+  Reset();
+
+  return kTRUE;
+}