Reordered ranom cones, now exclude areas close to the two leading jets
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetCluster.cxx
index bfdcd39..bbf947d 100644 (file)
@@ -21,7 +21,7 @@
 
  
 #include <TROOT.h>
-#include <TRandom.h>
+#include <TRandom3.h>
 #include <TSystem.h>
 #include <TInterpreter.h>
 #include <TChain.h>
@@ -55,7 +55,7 @@
 #include "AliJetKineReaderHeader.h"
 #include "AliGenCocktailEventHeader.h"
 #include "AliInputEventHandler.h"
-
+#include "AliAODJetEventBackground.h"
 
 #include "fastjet/PseudoJet.hh"
 #include "fastjet/ClusterSequenceArea.hh"
@@ -69,6 +69,7 @@ ClassImp(AliAnalysisTaskJetCluster)
 
 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
   delete fRef;
+  delete fRandom;
 }
 
 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
@@ -78,22 +79,32 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
   fUseAODTrackInput(kFALSE),
   fUseAODMCInput(kFALSE),
   fUseGlobalSelection(kFALSE),
+  fUseBackgroundCalc(kFALSE),
   fFilterMask(0),
   fTrackTypeRec(kTrackUndef),
   fTrackTypeGen(kTrackUndef),  
   fNSkipLeadingRan(0),
+  fNRandomCones(10),
   fAvgTrials(1),
   fExternalWeight(1),    
   fRecEtaWindow(0.5),
   fTrackPtCut(0.),                                                     
   fJetOutputMinPt(1),
+  fJetTriggerPtCut(0),
+  fCentCutUp(0),
+  fCentCutLo(0),
   fNonStdBranch(""),
+  fBackgroundBranch(""),
   fNonStdFile(""),
   fRparam(1.0), 
   fAlgorithm(fastjet::kt_algorithm),
   fStrategy(fastjet::Best),
   fRecombScheme(fastjet::BIpt_scheme),
   fAreaType(fastjet::active_area), 
+  fGhostArea(0.01),
+  fActiveAreaRepeats(1),
+  fGhostEtamax(1.5),
+  fRandom(0),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
@@ -115,6 +126,12 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
   fh1PtJetsRecInRan(0x0),
   fh1PtTracksGenIn(0x0),
   fh1Nch(0x0),
+  fh1CentralityPhySel(0x0), 
+  fh1Centrality(0x0), 
+  fh1CentralitySelect(0x0),
+  fh1ZPhySel(0x0), 
+  fh1Z(0x0), 
+  fh1ZSelect(0x0),
   fh2NRecJetsPt(0x0),
   fh2NRecTracksPt(0x0),
   fh2NConstPt(0x0),
@@ -144,7 +161,16 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
   fh2TracksLeadingJetPhiPtWRan(0x0),
   fHistList(0x0)  
 {
-
+  for(int i = 0;i<3;i++){
+    fh1BiARandomCones[i] = 0;
+    fh1BiARandomConesRan[i] = 0;
+  }
+  for(int i = 0;i<kMaxCent;i++){
+    fh2JetsLeadingPhiPtC[i] = 0;     
+    fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
+    fh2TracksLeadingJetPhiPtC[i] = 0;
+    fh2TracksLeadingJetPhiPtWC[i] = 0;
+  }
 }
 
 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
@@ -155,22 +181,32 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
   fUseAODTrackInput(kFALSE),
   fUseAODMCInput(kFALSE),
   fUseGlobalSelection(kFALSE),
+  fUseBackgroundCalc(kFALSE),
   fFilterMask(0),
   fTrackTypeRec(kTrackUndef),
   fTrackTypeGen(kTrackUndef),
   fNSkipLeadingRan(0),
+  fNRandomCones(5),
   fAvgTrials(1),
   fExternalWeight(1),    
   fRecEtaWindow(0.5),
   fTrackPtCut(0.),                                                     
   fJetOutputMinPt(1),
+  fJetTriggerPtCut(0),
+  fCentCutUp(0),
+  fCentCutLo(0),
   fNonStdBranch(""),
+  fBackgroundBranch(""),
   fNonStdFile(""),
   fRparam(1.0), 
   fAlgorithm(fastjet::kt_algorithm),
   fStrategy(fastjet::Best),
   fRecombScheme(fastjet::BIpt_scheme),
   fAreaType(fastjet::active_area), 
+  fGhostArea(0.01),
+  fActiveAreaRepeats(1),
+  fGhostEtamax(1.5),
+  fRandom(0),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
@@ -192,6 +228,12 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
   fh1PtJetsRecInRan(0x0),
   fh1PtTracksGenIn(0x0),
   fh1Nch(0x0),
+  fh1CentralityPhySel(0x0), 
+  fh1Centrality(0x0), 
+  fh1CentralitySelect(0x0),
+  fh1ZPhySel(0x0), 
+  fh1Z(0x0), 
+  fh1ZSelect(0x0),
   fh2NRecJetsPt(0x0),
   fh2NRecTracksPt(0x0),
   fh2NConstPt(0x0),
@@ -221,7 +263,17 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
   fh2TracksLeadingJetPhiPtWRan(0x0),
   fHistList(0x0)
 {
-  DefineOutput(1, TList::Class());  
+  for(int i = 0;i<3;i++){
+    fh1BiARandomCones[i] = 0;
+    fh1BiARandomConesRan[i] = 0;
+  }
+  for(int i = 0;i<kMaxCent;i++){
+    fh2JetsLeadingPhiPtC[i] = 0;     
+    fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
+    fh2TracksLeadingJetPhiPtC[i] = 0;
+    fh2TracksLeadingJetPhiPtWC[i] = 0;
+  }
+ DefineOutput(1, TList::Class());  
 }
 
 
@@ -242,7 +294,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
   // Create the output container
   //
 
-
+  fRandom = new TRandom3(0);
 
 
   // Connect the AOD
@@ -250,7 +302,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 
   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 
-
+  
 
   if(fNonStdBranch.Length()!=0)
     {
@@ -262,6 +314,36 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
       TClonesArray *tca = new TClonesArray("AliAODJet", 0);
       tca->SetName(fNonStdBranch.Data());
       AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
+
+   
+      TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
+      tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
+      AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
+
+
+      if(fUseBackgroundCalc){
+       if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
+         AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
+         evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
+         AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());  
+       }
+       // create the branch for the random cones with the same R 
+       TString cName = Form("%sRandomCone",fNonStdBranch.Data());
+       if(!AODEvent()->FindListObject(cName.Data())){
+         TClonesArray *tcaC = new TClonesArray("AliAODJet", 0);
+         tcaC->SetName(cName.Data());
+         AddAODBranch("TClonesArray",&tcaC,fNonStdFile.Data());
+       }
+       
+       // create the branch with the random for the random cones on the random event
+       cName = Form("%sRandomCone_random",fNonStdBranch.Data());
+       if(!AODEvent()->FindListObject(cName.Data())){
+         TClonesArray *tcaCran = new TClonesArray("AliAODJet", 0);
+         tcaCran->SetName(cName.Data());
+         AddAODBranch("TClonesArray",&tcaCran,fNonStdFile.Data());
+       }
+      }
+
       if(fNonStdFile.Length()!=0){
        // 
        // case that we have an AOD extension we need to fetch the jets from the extended output
@@ -287,8 +369,8 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
     }
 
 
-  OpenFile(1);
   if(!fHistList)fHistList = new TList();
+  fHistList->SetOwner();
 
   Bool_t oldStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
@@ -303,7 +385,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
       binLimitsPt[iPt] = 0.0;
     }
     else {// 1.0
-      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 0.5;
+      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
     }
   }
   
@@ -331,7 +413,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
     }
   }
 
-  const int nChMax = 100;
+  const int nChMax = 4000;
 
   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
@@ -364,6 +446,14 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 
+  fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
+  fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
+  fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
+
+  fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
+  fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
+  fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
+
   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
@@ -423,6 +513,18 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
                                       nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 
 
+  if(fUseBackgroundCalc){
+    for(int i = 0;i<3;i++){
+      fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
+      fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
+    }
+  }
+  for(int i = 0;i < kMaxCent;i++){
+    fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
+    fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
+    fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
+    fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
+  }
 
   const Int_t saveLevel = 3; // large save level more histos
   if(saveLevel>0){
@@ -444,6 +546,25 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
     fHistList->Add(fh1NConstLeadingRecRan);
     fHistList->Add(fh1PtJetsRecInRan);
     fHistList->Add(fh1Nch);
+    fHistList->Add(fh1Centrality);
+    fHistList->Add(fh1CentralitySelect);
+    fHistList->Add(fh1CentralityPhySel);
+    fHistList->Add(fh1Z);
+    fHistList->Add(fh1ZSelect);
+    fHistList->Add(fh1ZPhySel);
+    if(fUseBackgroundCalc){
+      for(int i = 0;i<3;i++){
+       fHistList->Add(fh1BiARandomCones[i]);
+       fHistList->Add(fh1BiARandomConesRan[i]);
+      }
+    }
+  for(int i = 0;i < kMaxCent;i++){
+    fHistList->Add(fh2JetsLeadingPhiPtC[i]);
+    fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
+    fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
+    fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
+  }
+
     fHistList->Add(fh2NRecJetsPt);
     fHistList->Add(fh2NRecTracksPt);
     fHistList->Add(fh2NConstPt);
@@ -506,15 +627,44 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
     return;
   }
 
+
+
   // handle and reset the output jet branch 
   // only need this once
   TClonesArray* jarray = 0;  
+  TClonesArray* jarrayran = 0;  
+  TClonesArray* rConeArray = 0;  
+  TClonesArray* rConeArrayRan = 0;  
+  AliAODJetEventBackground*  evBkg = 0;
   if(fNonStdBranch.Length()!=0) {
     if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
     if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
     if(jarray)jarray->Delete();    // this is our responsibility, clear before filling again
+    if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
+    if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
+    if(jarrayran)jarrayran->Delete();    // this is our responsibility, clear before filling again
+  
+    if(fUseBackgroundCalc){
+      if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
+      if(!evBkg)  evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
+      if(evBkg)evBkg->Reset(); 
+      
+      TString cName = Form("%sRandomCone",fNonStdBranch.Data());
+      if(AODEvent())rConeArray = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
+      if(!rConeArray)rConeArray =  (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
+      if(rConeArray)rConeArray->Delete();
+
+      cName = Form("%sRandomCone_random",fNonStdBranch.Data());
+      if(AODEvent())rConeArrayRan = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
+      if(!rConeArrayRan)rConeArrayRan =  (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
+      if(rConeArrayRan)rConeArrayRan->Delete();
+    }
   }
 
+  static AliAODJetEventBackground* externalBackground = 0;
+  if(!externalBackground&&fBackgroundBranch.Length()){
+    externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
+  }
   //
   // Execute analysis for current event
   //
@@ -542,25 +692,48 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
   Bool_t selectEvent =  false;
   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 
+  Float_t cent = 0;
+  Float_t zVtx  = 0;
+  Int_t cenClass = -1;
   if(fAOD){
     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
     TString vtxTitle(vtxAOD->GetTitle());
+    zVtx = vtxAOD->GetZ();
+
+    cent = fAOD->GetHeader()->GetCentrality();
+    if(cent<10)cenClass = 0;
+    else if(cent<30)cenClass = 1;
+    else if(cent<50)cenClass = 2;
+    else if(cent<80)cenClass = 3;
+    if(physicsSelection){
+      fh1CentralityPhySel->Fill(cent);
+      fh1ZPhySel->Fill(zVtx);
+    }
+
 
     if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
-       Float_t zvtx = vtxAOD->GetZ();
        Float_t yvtx = vtxAOD->GetY();
        Float_t xvtx = vtxAOD->GetX();
        Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
-       if(TMath::Abs(zvtx)<8.&&r2<1.){
-         if(physicsSelection)selectEvent = true;
+       if(TMath::Abs(zVtx)<20.&&r2<1.){ // apply vertex cut later on
+         if(physicsSelection){
+           selectEvent = true;
+         }
        }
     }
+    if(fCentCutUp>0){
+      if(cent<fCentCutLo||cent>fCentCutUp){
+       selectEvent = false;
+      }
+    }
+
   }
   if(!selectEvent){
     PostData(1, fHistList);
     return;
   }
-  
+  fh1Centrality->Fill(cent);  
+  fh1Z->Fill(zVtx);
   fh1Trials->Fill("#sum{ntrials}",1);
   
 
@@ -586,10 +759,14 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 
   vector<fastjet::PseudoJet> inputParticlesRec;
   vector<fastjet::PseudoJet> inputParticlesRecRan;
+  
+  // Generate the random cones
+  
+  AliAODJet vTmpRan(1,0,0,1);
   for(int i = 0; i < recParticles.GetEntries(); i++){
     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
-    // we talk total momentum here
+    // we take total momentum here
     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
     jInp.set_user_index(i);
     inputParticlesRec.push_back(jInp);
@@ -597,18 +774,20 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
     // the randomized input changes eta and phi, but keeps the p_T
     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
       Double_t pT = vp->Pt();
-      Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
-      Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
+      Double_t eta = 1.8 * fRandom->Rndm() - 0.9;
+      Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
       
-      Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));  
+      Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));  
       Double_t pZ = pT/TMath::Tan(theta);
 
       Double_t pX = pT * TMath::Cos(phi);
       Double_t pY = pT * TMath::Sin(phi);
       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
+
       jInpRan.set_user_index(i);
       inputParticlesRecRan.push_back(jInpRan);
+      vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
     }
 
     // fill the tref array, only needed when we write out jets
@@ -619,7 +798,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
       }
       fRef->Add(vp);
     }
-  }
+  }// recparticles
 
   if(inputParticlesRec.size()==0){
     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
@@ -628,10 +807,33 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
   }
   
   // run fast jet
+  // employ setters for these...
+
+  // now create the object that holds info about ghosts                        
+  if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
+    // reduce CPU time...
+    fGhostArea = 0.5; 
+    fActiveAreaRepeats = 0; 
+  }
+   
+ fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
+  fastjet::AreaType areaType =   fastjet::active_area;
+  fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
-  fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
+  fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
   
+  //range where to compute background
+  Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
+  phiMin = 0;
+  phiMax = 2*TMath::Pi();
+  rapMax = fGhostEtamax - fRparam;
+  rapMin = - fGhostEtamax + fRparam;
+  fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
+
+
   inclusiveJets = clustSeq.inclusive_jets();
   sortedJets    = sorted_by_pt(inclusiveJets);
  
@@ -643,6 +845,8 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
   Int_t nRec     = inclusiveJets.size();
   if(inclusiveJets.size()>0){
     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
+    Double_t area = clustSeq.area(sortedJets[0]);
+    leadingJet.SetEffArea(area,0);
     Float_t pt = leadingJet.Pt();
     Int_t nAodOutJets = 0;
     Int_t nAodOutTracks = 0;
@@ -664,8 +868,13 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
     Float_t phi = leadingJet.Phi();
     if(phi<0)phi+=TMath::Pi()*2.;    
     Float_t eta = leadingJet.Eta();
-    pt = leadingJet.Pt();
-    
+    Float_t pTback = 0;
+    if(externalBackground){
+      // carefull has to be filled in a task before
+      // todo, ReArrange to the botom
+      pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
+    }
+    pt = leadingJet.Pt() - pTback;
     // correlation of leading jet with tracks
     TIterator *recIter = recParticles.MakeIterator();
     recIter->Reset();
@@ -680,30 +889,157 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
+      if(cenClass>=0){
+       fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
+       fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
+      }
+
     }  
     
    
     
-   for(int j = 0; j < nRec;j++){
-     AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
-     aodOutJet = 0;
-     nAodOutTracks = 0;
-     Float_t tmpPt = tmpRec.Pt();  
-     fh1PtJetsRecIn->Fill(tmpPt);
-     // Fill Spectra with constituents
-     vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
-     fh1NConstRec->Fill(constituents.size());
-     fh2PtNch->Fill(nCh,tmpPt);
-     fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
-     fh2NConstPt->Fill(tmpPt,constituents.size());
-     // loop over constiutents and fill spectrum
-
-     // Add the jet information and the track references to the output AOD
-     
-     if(tmpPt>fJetOutputMinPt&&jarray){
-       aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
-     }
-
+    for(int j = 0; j < nRec;j++){
+      AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
+      aodOutJet = 0;
+      nAodOutTracks = 0;
+      Float_t tmpPt = tmpRec.Pt();  
+      Float_t tmpPtBack = 0;
+      if(externalBackground){
+       // carefull has to be filled in a task before
+       // todo, ReArrange to the botom
+       tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
+      }
+      tmpPt = tmpPt - tmpPtBack;
+      if(tmpPt<0)tmpPt = 0; // avoid negative weights...
+      
+      fh1PtJetsRecIn->Fill(tmpPt);
+      // Fill Spectra with constituents
+      vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
+      fh1NConstRec->Fill(constituents.size());
+      fh2PtNch->Fill(nCh,tmpPt);
+      fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
+      fh2NConstPt->Fill(tmpPt,constituents.size());
+      // loop over constiutents and fill spectrum
+      
+      // Add the jet information and the track references to the output AOD
+      
+      if(tmpPt>fJetOutputMinPt&&jarray){
+       aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
+       Double_t area1 = clustSeq.area(sortedJets[j]);
+       aodOutJet->SetEffArea(area1,0);
+      }
+            
+
+      if(fUseBackgroundCalc){       
+       
+       // create a random jet within the acceptance
+       Double_t etaMax = 0.9 - fRparam;
+       Int_t nCone = 0;
+       Int_t nConeRan = 0;
+       Double_t pTC = 1; // small number
+       for(int ir = 0;ir < fNRandomCones;ir++){
+         Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
+         Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
+         // massless jet
+         Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
+         Double_t pZC = pTC/TMath::Tan(thetaC);
+         Double_t pXC = pTC * TMath::Cos(phiC);
+         Double_t pYC = pTC * TMath::Sin(phiC);
+         Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
+         AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
+         bool skip = false;
+         for(int jj = 0; jj < TMath::Min(nRec,2);jj++){
+           AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
+           if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
+             skip = true;
+             break;
+           }
+         }
+         if(skip)continue;
+         tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
+         if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
+         if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
+       }// random cones
+
+
+       // loop over the reconstructed particles and add up the pT in the random cones
+       // maybe better to loop over randomized particles not in the real jets...
+       // but this by definition brings dow average energy in the whole  event
+       AliAODJet vTmpRanR(1,0,0,1);
+       for(int i = 0; i < recParticles.GetEntries(); i++){
+         AliVParticle *vp = (AliVParticle*)recParticles.At(i);
+         if(rConeArray){
+           for(int ir = 0;ir < fNRandomCones;ir++){
+             AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);  
+             if(jC&&jC->DeltaR(vp)<fRparam){
+               jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
+             }
+           }  
+         }// add up energy in cone
+      
+         // the randomized input changes eta and phi, but keeps the p_T
+         if(i>=fNSkipLeadingRan){// eventually skip the leading particles
+           Double_t pTR = vp->Pt();
+           Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
+           Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
+           
+           Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));  
+           Double_t pZR = pTR/TMath::Tan(thetaR);
+       
+           Double_t pXR = pTR * TMath::Cos(phiR);
+           Double_t pYR = pTR * TMath::Sin(phiR);
+           Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR); 
+           vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
+           if(rConeArrayRan){
+             for(int ir = 0;ir < fNRandomCones;ir++){
+               AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);  
+               if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
+                 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
+               }
+             }  
+           }
+         }
+       }// loop over recparticles
+    
+       Float_t jetArea = fRparam*fRparam*TMath::Pi();
+       for(int ir = 0;ir < fNRandomCones;ir++){
+         // rescale the momntum vectors for the random cones
+         if(!rConeArray)continue;
+         AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
+         if(rC){
+           Double_t etaC = rC->Eta();
+           Double_t phiC = rC->Phi();
+           // massless jet, unit vector
+           pTC = rC->ChargedBgEnergy();
+           if(pTC<=0)pTC = 0.1; // for almost empty events
+           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
+           Double_t pZC = pTC/TMath::Tan(thetaC);
+           Double_t pXC = pTC * TMath::Cos(phiC);
+           Double_t pYC = pTC * TMath::Sin(phiC);
+           Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
+           rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
+           rC->SetBgEnergy(0,0);
+           rC->SetEffArea(jetArea,0);
+         }
+         rC = (AliAODJet*)rConeArrayRan->At(ir);
+         // same wit random
+         if(rC){
+           Double_t etaC = rC->Eta();
+           Double_t phiC = rC->Phi();
+           // massless jet, unit vector
+           pTC = rC->ChargedBgEnergy();
+           if(pTC<=0)pTC = 0.1;// for almost empty events
+           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
+           Double_t pZC = pTC/TMath::Tan(thetaC);
+           Double_t pXC = pTC * TMath::Cos(phiC);
+           Double_t pYC = pTC * TMath::Sin(phiC);
+           Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
+           rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
+           rC->SetBgEnergy(0,0);
+           rC->SetEffArea(jetArea,0);
+         }
+       }
+      }// if(fUseBackgroundCalc
      
      for(unsigned int ic = 0; ic < constituents.size();ic++){
        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
@@ -713,19 +1049,14 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        }
        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
      }
-
-
-
-     
-     
-
-
      
      // correlation
      Float_t tmpPhi =  tmpRec.Phi();
      Float_t tmpEta =  tmpRec.Eta();
      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
      
+
+     
      if(j==0){
        fh1PtJetsLeadingRecIn->Fill(tmpPt);
        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
@@ -742,18 +1073,54 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
      Float_t dEta = eta - tmpRec.Eta();
      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
+     if(cenClass>=0){
+       fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
+       fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
+     }
      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
-   }
+    }
    delete recIter;
+  
+   //background estimates:all bckg jets(0) & wo the 2 hardest(1)
+
+   if(evBkg){
+     vector<fastjet::PseudoJet> jets2=sortedJets;
+     if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); 
+     Double_t bkg1=0;
+     Double_t sigma1=0.;
+     Double_t meanarea1=0.;
+     Double_t bkg2=0;
+     Double_t sigma2=0.;
+     Double_t meanarea2=0.;
+
+     clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
+     evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+
+     //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));    
+     //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));    
+     
+     clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
+     evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
+     //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));    
+     //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));    
+
+   }
   }
+   
+
+  
+  
  
- // fill track information
- Int_t nTrackOver = recParticles.GetSize();
+  // fill track information
+  Int_t nTrackOver = recParticles.GetSize();
   // do the same for tracks and jets
- if(nTrackOver>0){
+
+  if(nTrackOver>0){
    TIterator *recIter = recParticles.MakeIterator();
    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
    Float_t pt = tmpRec->Pt();
+
    //    Printf("Leading track p_t %3.3E",pt);
    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
@@ -801,7 +1168,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 
  // find the random jets
  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
- fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
+ fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
   
  inclusiveJetsRan = clustSeqRan.inclusive_jets();
  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
@@ -850,8 +1217,8 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
       }  
 
-
-
+    Int_t nAodOutJetsRan = 0;
+     AliAODJet *aodOutJetRan = 0;
     for(int j = 0; j < nRecRan;j++){
       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
       Float_t tmpPt = tmpRec.Pt();
@@ -862,6 +1229,28 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
       fh2NConstPtRan->Fill(tmpPt,constituents.size());
       fh2PtNchRan->Fill(nCh,tmpPt);
       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
+
+
+     if(tmpPt>fJetOutputMinPt&&jarrayran){
+       aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
+       Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
+       
+       aodOutJetRan->SetEffArea(arearan,0);    }
+
+
+
+     
+     for(unsigned int ic = 0; ic < constituents.size();ic++){
+       // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
+       // fh1PtJetConstRec->Fill(part->Pt());
+       if(aodOutJetRan){
+        aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
+       }
+     }
+      
+
+
+
       // correlation
       Float_t tmpPhi =  tmpRec.Phi();
       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
@@ -873,8 +1262,58 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        continue;
       }
     }  
+
+     
+    if(evBkg){
+     Double_t bkg3=0.;
+     Double_t sigma3=0.;
+     Double_t meanarea3=0.;
+     clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
+     evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
+     //     float areaRandomCone = rRandomCone2 *TMath::Pi();         
+     /*
+     fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));    
+     fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));    
+     */
+    }
+
+
+
+ }
+
+
+ // do the event selection if activated
+ if(fJetTriggerPtCut>0){
+   bool select = false;
+   Float_t minPt = fJetTriggerPtCut;
+   /*
+   // hard coded for now ...
+   // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
+   if(cent<10)minPt = 50;
+   else if(cent<30)minPt = 42;
+   else if(cent<50)minPt = 28;
+   else if(cent<80)minPt = 18;
+   */
+   float rho = 0;
+   if(externalBackground)rho = externalBackground->GetBackground(2);
+   if(jarray){
+     for(int i = 0;i < jarray->GetEntriesFast();i++){
+       AliAODJet *jet = (AliAODJet*)jarray->At(i);
+       Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
+       if(ptSub>=minPt){
+        select = true;
+        break;
+       }
+     }
+   }   
+
+   if(select){
+     static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+     fh1CentralitySelect->Fill(cent);
+     fh1ZSelect->Fill(zVtx);
+     aodH->SetFillAOD(kTRUE);
+   }
  }
 
  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
  PostData(1, fHistList);
@@ -904,7 +1343,7 @@ Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
       AliAODTrack *tr = aod->GetTrack(it);
       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
       if(TMath::Abs(tr->Eta())>0.9)continue;
-      //      if(tr->Pt()<0.3)continue;
+      if(tr->Pt()<fTrackPtCut)continue;
       list->Add(tr);
       iCount++;
     }
@@ -917,11 +1356,13 @@ Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
       if(!mcEvent->IsPhysicalPrimary(it))continue;
       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
       if(type == kTrackKineAll){
+       if(part->Pt()<fTrackPtCut)continue;
        list->Add(part);
        iCount++;
       }
       else if(type == kTrackKineCharged){
        if(part->Particle()->GetPDG()->Charge()==0)continue;
+       if(part->Pt()<fTrackPtCut)continue;
        list->Add(part);
        iCount++;
       }
@@ -938,11 +1379,13 @@ Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
       if(!part->IsPhysicalPrimary())continue;
       if(type == kTrackAODMCAll){
+       if(part->Pt()<fTrackPtCut)continue;
        list->Add(part);
        iCount++;
       }
       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
        if(part->Charge()==0)continue;
+       if(part->Pt()<fTrackPtCut)continue;
        if(kTrackAODMCCharged){
          list->Add(part);
        }
@@ -956,7 +1399,6 @@ Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
   }// AODMCparticle
   list->Sort();
   return iCount;
-
 }
 
 /*