new UE Task from Sara
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 Jun 2010 12:37:32 +0000 (12:37 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 Jun 2010 12:37:32 +0000 (12:37 +0000)
PWG4/JetTasks/AliAnalysisTaskUE.cxx
PWG4/JetTasks/AliAnalysisTaskUE.h

index b42fec1..5e07c4b 100644 (file)
 
 /* $Id:$ */
 
+//class AliAODInputHandler;
+//class AliAODHandler;
+//class AliLog;
+//class TROOT;
+
 #include <TROOT.h>
-#include <TSystem.h>
 #include <TChain.h>
 #include <TFile.h>
 #include <TList.h>
-#include <TH1I.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TProfile.h>
-#include <TCanvas.h>
-#include <TVector3.h>
-#include <TLorentzVector.h>
 #include <TMath.h>
 #include <TTree.h>
-#include <TBranch.h>
-#include <TRandom.h>
+#include <TVector3.h>
+#include <TH1F.h>
 
+#include "AliAnalyseUE.h"
 #include "AliAnalysisTaskUE.h"
+#include "AliHistogramsUE.h"
+
 #include "AliAnalysisManager.h"
 #include "AliMCEventHandler.h"
-#include "AliMCEvent.h"
-#include "AliAODEvent.h"
-#include "AliAODInputHandler.h"
-#include "AliAODHandler.h"
-#include "AliStack.h"
-#include "AliAODJet.h"
-#include "AliAODTrack.h"
-#include "AliAODMCParticle.h"
-#include "AliKFVertex.h"
 
-#include "AliGenPythiaEventHeader.h"
 #include "AliAnalysisHelperJetTasks.h"
-#include "AliInputEventHandler.h"
-#include "AliStack.h"
+#include "AliAODHandler.h"
+#include "AliAODInputHandler.h"
+#include "AliGenPythiaEventHeader.h"
 #include "AliLog.h"
+#include "AliInputEventHandler.h"
 
-//__________________________________________________________________
+////////////////////////////////////////////////////////////////////////
+//
 // Analysis class for Underlying Event studies
 //
 // Look for correlations on the tranverse regions to 
 //    Arian.Abrahantes.Quintana@cern.ch 
 //    Ernesto.Lopez.Torres@cern.ch
 // 
+////////////////////////////////////////////////////////////////////////
 
 
 ClassImp( AliAnalysisTaskUE)
 
-////////////////////////////////////////////////////////////////////////
-
+// Define global pointer
+AliAnalysisTaskUE* AliAnalysisTaskUE::fgTaskUE=NULL;
 
 //____________________________________________________________________
 AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
 AliAnalysisTask(name,""),
-fTrigger(0),
-fDebug(0),
-fDeltaAOD(kFALSE),
-fDeltaAODBranch(""),
-fAODBranch("jets"),
-fArrayJets(0x0),           
+fAnaUE(0x0),
 fAOD(0x0),            
-fAODjets(0x0),
+fAODBranch("jets"),
+fDebug(0),
 fListOfHistos(0x0),  
 fBinsPtInHist(30),     
-fMinJetPtInHist(0.),
-fMaxJetPtInHist(300.), 
 fIsNorm2Area(kTRUE),
-fUseMCParticleBranch(kFALSE),
+fMaxJetPtInHist(300.), 
+fMinJetPtInHist(0.),
 fConstrainDistance(kTRUE),
 fMinDistance(0.2),
 fSimulateChJetPt(kFALSE),
 fUseAliStack(kTRUE),
+fUseMCParticleBranch(kFALSE),
 fnTracksVertex(3),  // QA tracks pointing to principal vertex (= 3 default) 
 fZVertex(5.),
 fAnaType(1),         
-fRegionType(1),
-fConeRadius(0.7),
 fConePosition(1),
-fAreaReg(1.5393), // Pi*0.7*0.7
-fUseChPartJet(kFALSE),
-fUseChPartMaxPt(kFALSE),
+fConeRadius(0.7),
+fFilterBit(0xFF),
+fJetsOnFly(kFALSE),
+fRegionType(1),
 fUseChargeHadrons(kFALSE),
-fUseSingleCharge(kFALSE),
+fUseChPartJet(kFALSE),
 fUsePositiveCharge(kTRUE),
+fUseSingleCharge(kFALSE),
 fOrdering(1),
-fFilterBit(0xFF),
-fJetsOnFly(kFALSE),
 fChJetPtMin(5.0),
 fJet1EtaCut(0.2),
 fJet2DeltaPhiCut(2.616),    // 150 degrees
 fJet2RatioPtCut(0.8),
 fJet3PtCut(15.),
-fTrackPtCut(0.),
 fTrackEtaCut(0.9),
-fAvgTrials(1),
-fhNJets(0x0),
-fhEleadingPt(0x0),
-fhMinRegPtDist(0x0),
-fhRegionMultMin(0x0),
-fhMinRegAvePt(0x0), 
-fhMinRegSumPt(0x0),            
-fhMinRegMaxPtPart(0x0),
-fhMinRegSumPtvsMult(0x0),
-fhdNdEtaPhiDist(0x0),        
-fhFullRegPartPtDistVsEt(0x0), 
-fhTransRegPartPtDistVsEt(0x0),
-fhRegionSumPtMaxVsEt(0x0),
-fhRegionMultMax(0x0),         
-fhRegionMultMaxVsEt(0x0),     
-fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),         
-fhRegionMultMinVsEt(0x0),
-fhRegionAveSumPtVsEt(0x0),
-fhRegionDiffSumPtVsEt(0x0),
-fhRegionAvePartPtMaxVsEt(0x0),
-fhRegionAvePartPtMinVsEt(0x0),
-fhRegionMaxPartPtMaxVsEt(0x0),
-//fhRegForwardSumPtVsEt(0x0),
-//fhRegForwardMultVsEt(0x0),
-//fhRegBackwardSumPtVsEt(0x0),
-//fhRegBackwardMultVsEt(0x0),
-fhRegForwardMult(0x0),
-fhRegForwardSumPtvsMult(0x0),
-fhRegBackwardMult(0x0),
-fhRegBackwardSumPtvsMult(0x0),
-fhRegForwardPartPtDistVsEt(0x0),
-fhRegBackwardPartPtDistVsEt(0x0),
-fhRegTransMult(0x0),
-fhRegTransSumPtVsMult(0x0),
-fhMinRegSumPtJetPtBin(0x0),
-fhMaxRegSumPtJetPtBin(0x0),
-fhVertexMult(0x0),
-fh1Xsec(0x0),
-fh1Trials(0x0),
-fSettingsTree(0x0)//,   fhValidRegion(0x0)
+fTrackPtCut(0.),
+fAvgTrials(1)
 {
   // Default constructor
   // Define input and output slots here
@@ -170,6 +119,56 @@ fSettingsTree(0x0)//,   fhValidRegion(0x0)
 
 }
 
+//____________________________________________________________________
+AliAnalysisTaskUE:: AliAnalysisTaskUE(const AliAnalysisTaskUE & original):
+AliAnalysisTask(),
+fAnaUE(original.fAnaUE),
+fAOD(original.fAOD),            
+fAODBranch(original.fAODBranch),
+fDebug(original.fDebug),
+fListOfHistos(original.fListOfHistos),  
+fBinsPtInHist(original.fBinsPtInHist),     
+fIsNorm2Area(original.fIsNorm2Area),
+fMaxJetPtInHist(original.fMaxJetPtInHist), 
+fMinJetPtInHist(original.fMinJetPtInHist),
+fConstrainDistance(original.fConstrainDistance),
+fMinDistance(original.fMinDistance),
+fSimulateChJetPt(original.fSimulateChJetPt),
+fUseAliStack(original.fUseAliStack),
+fUseMCParticleBranch(original.fUseMCParticleBranch),
+fnTracksVertex(original.fnTracksVertex),  // QA tracks pointing to principal vertex (= 3 default) 
+fZVertex(original.fZVertex),
+fAnaType(original.fAnaType),         
+fConePosition(original.fConePosition),
+fConeRadius(original.fConeRadius),
+fFilterBit(original.fFilterBit),
+fJetsOnFly(original.fJetsOnFly),
+fRegionType(original.fRegionType),
+fUseChargeHadrons(original.fUseChargeHadrons),
+fUseChPartJet(original.fUseChPartJet),
+fUsePositiveCharge(original.fUsePositiveCharge),
+fUseSingleCharge(original.fUseSingleCharge),
+fOrdering(original.fOrdering),
+fChJetPtMin(original.fChJetPtMin),
+fJet1EtaCut(original.fJet1EtaCut),
+fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),    // 150 degrees
+fJet2RatioPtCut(original.fJet2RatioPtCut),
+fJet3PtCut(original.fJet3PtCut),
+fTrackEtaCut(original.fTrackEtaCut),
+fTrackPtCut(original.fTrackPtCut),
+fAvgTrials(original.fAvgTrials)
+{
+  // Copy constructor
+}
+
+
+//______________________________________________________________
+AliAnalysisTaskUE & AliAnalysisTaskUE::operator = (const AliAnalysisTaskUE & /*source*/)
+{
+  // assignment operator
+  return *this;
+}
+
 //______________________________________________________________
 Bool_t AliAnalysisTaskUE::Notify()
 {
@@ -180,24 +179,22 @@ Bool_t AliAnalysisTaskUE::Notify()
   fAvgTrials = 1;
   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
   Float_t xsection = 0;
-  Float_t ftrials  = 1;
+  Float_t trials  = 1;
   if(tree){
-    TFile *curfile = tree->GetCurrentFile();
-    if (!curfile) {
-      Error("Notify","No current file");
-      return kFALSE;
-    }
-    if(!fh1Xsec||!fh1Trials){
-      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
-      return kFALSE;
-    }
-    AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials); 
-    fh1Xsec->Fill("<#sigma>",xsection);
+       TFile *curfile = tree->GetCurrentFile();
+       if (!curfile) {
+               Error("Notify","No current file");
+               return kFALSE;
+               }
+       
+       AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials); 
+        
+       fAnaUE->FillXsec("<#sigma>",xsection);
 
-   // construct average trials
-   Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
-   if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
-  }
+       // construct average trials
+       Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
+       if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
+       }
   
   return kTRUE;
 }
@@ -212,7 +209,7 @@ void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
   // or to the OutputEventHandler ( AOD is create by a previus task in the train )
   // we need to check where it is and get the pointer to AODEvent in the right way
   
-  // Delta AODs are also implemented
+  // Delta AODs are also accepted
   
  
   if (fDebug > 1) AliInfo("ConnectInputData() ");
@@ -220,47 +217,44 @@ void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
   
   if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD
-    fAOD = ((AliAODInputHandler*)handler)->GetEvent();
-    if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler");
-    // Case when jets are reconstructed on the fly from AOD tracks
-    // (the Jet Finder is using the AliJetAODReader) of InputEventHandler
-    // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with
-    // different parameters to default ones stored in the AOD or to use a different algorithm
-    if( fJetsOnFly ) {
-      handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
-      if( handler && handler->InheritsFrom("AliAODHandler") ) {
-        fAODjets = ((AliAODHandler*)handler)->GetAOD();
-        if (fDebug > 1) AliInfo(" ==== Jets from AliAODHandler (on the fly)");
-      }
-    } else {
-      fAODjets = fAOD;
-      if (fDebug > 1) AliInfo(" ==== Jets from AliAODInputHandler");
-    }
-  } else {  //output AOD
-    handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
-    if( handler && handler->InheritsFrom("AliAODHandler") ) {
-      fAOD = ((AliAODHandler*)handler)->GetAOD();
-      fAODjets = fAOD;
-      if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
-    } else {
-      AliFatal("I can't get any AOD Event Handler");
-      return;
-    }
-  }
+       fAOD = ((AliAODInputHandler*)handler)->GetEvent();
+       if(!fJetsOnFly){
+               if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
+               }else{
+               AliInfo(" ==== Tracks from AliAODInputHandler / Jets on-the-fly");
+               }
+       } else {  //output AOD
+       handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
+       if( handler && handler->InheritsFrom("AliAODHandler") ) {
+               fAOD = ((AliAODHandler*)handler)->GetAOD();
+               if (!fJetsOnFly){
+                       if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
+                       } else {
+                       if (fDebug > 1) AliInfo(" ==== Tracks from AliAODHandler / Jets on-the-fly");
+                       }
+               }else {
+               AliFatal("I can't get any AOD Event Handler");
+               return;
+               }
+       }       
+
+   fAnaUE->Initialize( *this );
 }
 
 //____________________________________________________________________
 void  AliAnalysisTaskUE::CreateOutputObjects()
 {
   // Create the output container
-  //
+  
   if (fDebug > 1) AliInfo("CreateOutPutData()");
-  //
-  //  Histograms
-
-  // OpenFile(0);
-  CreateHistos();
-  fListOfHistos->SetOwner(kTRUE);  
+   
+  //Initialize AliAnalysisUE, a class implementing the main analysis algorithms
+  fAnaUE = new AliAnalyseUE();
+  
+  //Initialize output histograms
+  fAnaUE->CreateHistograms(fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
+  fListOfHistos = (TList*)fAnaUE->GetHistograms()->Clone();
+  fListOfHistos->SetOwner(kTRUE); 
   PostData(0,fListOfHistos);
 
 }
@@ -268,335 +262,76 @@ void  AliAnalysisTaskUE::CreateOutputObjects()
 //____________________________________________________________________
 void  AliAnalysisTaskUE::Exec(Option_t */*option*/)
 {
-  //Trigger selection ************************************************
-  /*AliInputEventHandler* inputHandler = (AliInputEventHandler*)
+  // Trigger selection ************************************************
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*)
          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
-  if (inputHandler->IsEventSelected()) {
-    if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
-  } else {
-    if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
-    return;
-  }
-   */                                     
-  //Event selection (vertex) *****************************************
-
-  Int_t nVertex = fAOD->GetNumberOfVertices();
-  if( nVertex > 0 ) { // Only one vertex (reject pileup??)
-     AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
-     TString tpcvertex("TPCVertex");
-     if (vertex->GetName() == tpcvertex){
-       if (fDebug > 1) AliInfo(Form("Primary vertex selection: %s event REJECTED ...",vertex->GetName()));
-       return;
-        }
-     Int_t nTracksPrim = vertex->GetNContributors();
-     Double_t zVertex = vertex->GetZ();
-     if (fDebug > 1) AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
-     // Select a quality vertex by number of tracks?
-     if( nTracksPrim < fnTracksVertex || TMath::Abs(zVertex) > fZVertex ) {
-        if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
-        return;
-     }
-     if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
-  } else {
-     if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
-     return;
-  } 
-  /*AliKFVertex primVtx(*(fAOD->GetPrimaryVertex()));
-  Int_t nTracksPrim=primVtx.GetNContributors();
-  if (fDebug > 1) AliInfo(Form(" Primary-vertex Selection: %d",nTracksPrim));
-  if(!nTracksPrim){
-    if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
-    return;
-  }
-  if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED ...");
-  */
-  // Execute analysis for current event
-  //
+  if( !inputHandler->InheritsFrom("AliAODInputHandler") ) { // input AOD
+       if (inputHandler->IsEventSelected()) {
+               if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
+       } else {
+               if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
+               return;
+               }
+        }                                
+  // Event selection (vertex) *****************************************
+   
+  if(!fAnaUE->VertexSelection(fAOD,fnTracksVertex,fZVertex)) return;
+  
+  // Execute analysis for current event ******************************
+  
   if ( fDebug > 3 ) AliInfo( " Processing event..." );
+
   // fetch the pythia header info and get the trials
   AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
   Float_t nTrials = 1;
   if (mcHandler) {
-       AliMCEvent* mcEvent = mcHandler->MCEvent();
-       if (mcEvent) {
-         AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
-         if(pythiaGenHeader){
-           nTrials = pythiaGenHeader->Trials();
-         }
-      }
-    }
-  fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
+       AliMCEvent* mcEvent = mcHandler->MCEvent();
+       if (mcEvent) {
+               AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
+               if(pythiaGenHeader){
+                       nTrials = pythiaGenHeader->Trials();
+                       }
+               }
+       }
+  fAnaUE->FillTrials("#sum{ntrials}",fAvgTrials);
+  
+  //analyse the event
   AnalyseUE();
-
-  // Post the data
-  PostData(0, fListOfHistos);
+  fListOfHistos = (TList*)fAnaUE->GetHistograms()->Clone();
+  PostData(0,fListOfHistos);
 }
 
 //____________________________________________________________________
 void  AliAnalysisTaskUE::AnalyseUE()
 {
-   Double_t const kMyTolerance = 0.0000001; 
-  //
-  // Look for correlations on the tranverse regions to 
-  // the leading charged jet
-  // 
-  
-  // ------------------------------------------------
-  // Find Leading Jets 1,2,3 
-  // (could be skipped if Jets are sort by Pt...)
-  Double_t maxPtJet1 = 0.; 
-  Int_t    index1 = -1;
-  Double_t maxPtJet2 = 0.;   // jet 2 need for back to back inclusive
-  Int_t    index2 = -1;
-  Double_t maxPtJet3 = 0.;   // jet 3 need for back to back exclusive
-  Int_t    index3 = -1;
-  TVector3 jetVect[3];
-  Int_t nJets = 0;
+
   
+   // Get jets and order by pT
+   TVector3 jetVect[3];
+   *jetVect = fAnaUE->GetOrderedClusters(fAODBranch, fUseChPartJet, fChJetPtMin);
   
-  if( !fUseChPartJet && fAnaType != 4 ) {
-    
-    // Use AOD Jets
-    if(fDeltaAOD){
-      if (fDebug > 1) AliInfo(" ==== Jets From  Delta-AODs !");
-      if (fDebug > 1) AliInfo(Form(" ====  Reading Branch: %s  ",fDeltaAODBranch.Data()));
-          fArrayJets = (TClonesArray*)fAODjets->GetList()->FindObject(fDeltaAODBranch.Data());
-    } else {
-      if (fDebug > 1) AliInfo(" ==== Read Standard-AODs  !");
-      if (fDebug > 1) AliInfo(Form(" ====  Reading Branch: %s  ",fAODBranch.Data()));
-      fArrayJets = (TClonesArray*)fAODjets->GetList()->FindObject(fAODBranch.Data());
-    }
-    if (!fArrayJets){
-       AliFatal(" No jet-array! ");
-       return;
-    }
-    nJets=fArrayJets->GetEntries();
-    //printf("AOD %d jets \n", nJets);
+   if( jetVect[0].Pt() < 0. ) {
+       if( fDebug > 1 ) AliInfo("\n   Skipping Event, not jet found...");
+       return;
+       } else {
+               if (fDebug >1 ) AliInfo(Form("\n   Pt Leading Jet = %6.1f eta=%5.3f ",  jetVect[0].Pt(), jetVect[0].Eta() ));
+               }
 
-    for( Int_t i=0; i<nJets; ++i ) {
-      AliAODJet* jet = (AliAODJet*)fArrayJets->At(i);
-      Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!!
-      if( jetPt > maxPtJet1 ) {
-            maxPtJet3 = maxPtJet2; index3 = index2;
-            maxPtJet2 = maxPtJet1; index2 = index1;
-            maxPtJet1 = jetPt; index1 = i;
-      } else if( jetPt > maxPtJet2 ) {
-            maxPtJet3 = maxPtJet2; index3 = index2;
-            maxPtJet2 = jetPt; index2 = i;
-      } else if( jetPt > maxPtJet3 ) {
-            maxPtJet3 = jetPt; index3 = i;
-      }
-    }
-
-    if( index1 != -1 ) {
-      AliAODJet *jet =(AliAODJet*) fArrayJets->At(index1);
-      if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
-    }
-    if( index2 != -1 ) {
-      AliAODJet* jet= (AliAODJet*) fArrayJets->At(index2);
-      if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
-    }
-    if( index3 != -1 ) {
-       AliAODJet* jet = (AliAODJet*) fArrayJets->At(index3);
-      if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
-    }
-    
-  }
-
-  if( fUseChPartJet )  {
-
-    // Use "Charged Particle Jets"
-    TObjArray* jets = FindChargedParticleJets();
-    if( jets ) {
-      nJets = jets->GetEntriesFast();
-      if( nJets > 0 ) {
-            index1 = 0;
-            AliAODJet* jet = (AliAODJet*)jets->At(0);
-            maxPtJet1 = jet->Pt();
-            jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
-      }
-      if( nJets > 1 ) {
-            index2 = 1;
-            AliAODJet* jet = (AliAODJet*)jets->At(1);
-        maxPtJet2 = jet->Pt();
-            jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
-      }
-      if( nJets > 2 ) {
-            index3 = 2;
-            AliAODJet* jet = (AliAODJet*)jets->At(2);
-        maxPtJet3 = jet->Pt();
-            jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
-      }
-      
-      jets->Delete();
-      delete jets;
-      if (fDebug > 1) AliInfo(" ==== Jets Created in OUR class  !: CDF algorithm ");
-    }
-  }
-
-
-  if( fAnaType == 4 )  {
-
-    // Use "Max Pt Charge Particle"
-    TObjArray* tracks = SortChargedParticles();
-    if( tracks ) {
-      nJets = tracks->GetEntriesFast();
-      if( nJets > 0 ) {
-        index1 = 0;
-        AliAODTrack* jet = (AliAODTrack*)tracks->At(0);
-        maxPtJet1 = jet->Pt();
-        jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
-      }
-      tracks->Clear();
-      delete tracks; 
-    }
-  }
+   // Select events according to analysis type
+   if ( ! (fAnaUE->AnaTypeSelection(jetVect))) return;
 
+  // Find max and min regions with real tracks
+  if (!fUseMCParticleBranch){
+       // Primary vertex distribution
+       AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
+       fAnaUE->FillVertex(vertex->GetNContributors());
+       
+       // Fill leading "jet" histogram
+       fAnaUE->FillLeadingJet(jetVect[0].Pt());
 
-  fhNJets->Fill(nJets);
-  
-  if( fDebug > 1 ) {
-    if( index1 < 0 ) {
-      AliInfo("\n   Skipping Event, not jet found...");
-      return;
-    } else
-      AliInfo(Form("\n   Pt Leading Jet = %6.1f eta=%5.3f ",  maxPtJet1, jetVect[0].Eta() ));
-  }
-  
-  
-  // ----------------------------------------------
-  // Cut events by jets topology
-  // fAnaType:
-  //     1 = inclusive,
-  //         - Jet1 |eta| < fJet1EtaCut
-  //     2 = back to back inclusive
-  //         - fulfill case 1
-  //         - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut
-  //         - Jet2.Pt/Jet1Pt > fJet2RatioPtCut
-  //     3 = back to back exclusive
-  //         - fulfill case 2
-  //         - Jet3.Pt < fJet3PtCut
-  
-  if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) {
-    if( fDebug > 1 ) AliInfo("\n   Skipping Event...Jet1 |eta| > fJet1EtaCut");
-    return;
-  }
-  // back to back inclusive
-  if( fAnaType > 1 && fAnaType < 4 && index2 == -1 ) {
-    if( fDebug > 1 ) AliInfo("\n   Skipping Event... no second Jet found");
-    return;
-  }
-  if( fAnaType > 1 && fAnaType < 4 && index2 > -1 ) {
-    if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
-       maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
-      if( fDebug > 1 ) AliInfo("\n   Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
-      return;
-    }
-  }
-  // back to back exclusive
-  if( fAnaType > 2 && fAnaType < 4 && index3 > -1 ) {
-    if( maxPtJet3 > fJet3PtCut ) {
-      if( fDebug > 1 ) AliInfo("\n   Skipping Event... Jet3.Pt > fJet3PtCut ");
-      return;
-    }
-  }
-  
-  //fhEleadingPt->Fill( maxPtJet1 );
+       fAnaUE->FindMaxMinRegions( jetVect,  fConePosition );
 
-  // ----------------------------------------------
-  // Find max and min regions
-  Double_t sumPtRegionPosit = 0.;
-  Double_t sumPtRegionNegat = 0.;
-  Double_t sumPtRegionForward = 0;
-  Double_t sumPtRegionBackward = 0;
-  Double_t maxPartPtRegion  = 0.;
-  Int_t    nTrackRegionPosit = 0;
-  Int_t    nTrackRegionNegat = 0;
-  Int_t    nTrackRegionForward = 0;
-  Int_t    nTrackRegionBackward = 0;
-  static Double_t const  kPI     = TMath::Pi();
-  static Double_t const  kTWOPI  = 2.*kPI;
-  static Double_t const  k270rad = 270.*kPI/180.;
-  static Double_t const  k120rad = 120.*kPI/180.;
-  
-  //Area for Normalization Purpose at Display histos
-  // Forward and backward
-  Double_t normArea = 1.;
-  // Transverse
-  if (fIsNorm2Area) {
-    SetRegionArea(jetVect);
-    normArea =  2.*fTrackEtaCut*k120rad ;
-  } else fAreaReg = 1.;
-  
-  if (!fUseMCParticleBranch){
-    AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
-    fhVertexMult->Fill( vertex->GetNContributors() );
-    
-    fhEleadingPt->Fill( maxPtJet1 );
-    Int_t nTracks = fAOD->GetNTracks();
-    if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
-    
-    for (Int_t ipart=0; ipart<nTracks; ++ipart) {
-      
-      AliAODTrack* part = fAOD->GetTrack( ipart );
-      if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
-      if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
-      if (!part->IsPrimaryCandidate()) continue; // reject whatever is not linked to collision point
-      // PID Selection: Reject everything but hadrons
-      Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion || 
-                        part->GetMostProbablePID()==AliAODTrack::kKaon || 
-                        part->GetMostProbablePID()==AliAODTrack::kProton;
-      if ( fUseChargeHadrons && !isHadron ) continue;
-      
-      if ( !part->Charge() ) continue; //Only charged
-      if ( fUseSingleCharge ) { // Charge selection
-        if ( fUsePositiveCharge && part->Charge() < 0.) continue; // keep Positives
-        if ( !fUsePositiveCharge && part->Charge() > 0.) continue; // keep Negatives
-      }
-      
-      if ( part->Pt() < fTrackPtCut ) continue;
-      if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
-      
-      TVector3 partVect(part->Px(), part->Py(), part->Pz());
-      Bool_t isFlagPart = kTRUE;
-      Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
-      if( deltaPhi > kTWOPI )  deltaPhi-= kTWOPI;
-      if (fAnaType != 4 ) fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
-      else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
-         fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
-         isFlagPart = kFALSE;
-      }
-      
-      fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
-      
-      Int_t region = IsTrackInsideRegion( jetVect, &partVect );  
-      
-      if (region == 1) {
-        if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
-        sumPtRegionPosit += part->Pt();
-        nTrackRegionPosit++;
-        fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
-      }
-      if (region == -1) {
-        if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
-        sumPtRegionNegat += part->Pt();
-        nTrackRegionNegat++;
-        fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
-      }
-      if (region == 2){ //forward
-        sumPtRegionForward += part->Pt();
-        nTrackRegionForward++;
-        fhRegForwardPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
-      }
-      if (region == -2){ //backward
-        sumPtRegionBackward += part->Pt();
-        nTrackRegionBackward++;
-        fhRegBackwardPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
-      }
-    }//end loop AOD tracks
-  }
-  else {
+  }else { 
     
     // this is the part we only use when we have MC information
     // More than a test for values of it also resumes the reconstruction efficiency of jets
@@ -620,1012 +355,57 @@ void  AliAnalysisTaskUE::AnalyseUE()
       return;
     }
     
-    //Get Jets from MC header
-    Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
-    AliAODJet pythiaGenJets[4];
-    TVector3 jetVectnew[4];
-    Int_t iCount = 0;
-    for(int ip = 0;ip < nPythiaGenJets;++ip){
-      if (iCount>3) break;
-      Float_t p[4];
-      pythiaGenHeader->TriggerJet(ip,p);
-      TVector3 tempVect(p[0],p[1],p[2]);
-      if ( TMath::Abs(tempVect.Eta())>fJet1EtaCut ) continue;
-      pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
-      jetVectnew[iCount].SetXYZ(pythiaGenJets[iCount].Px(), pythiaGenJets[iCount].Py(), pythiaGenJets[iCount].Pz());
-      iCount++;
-    }
-    
-    if (!iCount) return;// no jet in eta acceptance
-    
-    //Search the index of the nearest MC jet to the leading jet reconstructed from the input data
-    Int_t index = 0;
-    if (fConstrainDistance){
-      Float_t deltaR = 0.;
-      Float_t dRTemp = 0.;
-      for (Int_t i=0; i<iCount; i++){
-         if (!i) {
-         dRTemp = jetVectnew[i].DeltaR(jetVect[0]);
-         index = i;
-         }
-         deltaR = jetVectnew[i].DeltaR(jetVect[0]);
-         if (deltaR < dRTemp){
-         index = i;
-         dRTemp = deltaR;
-         }
-      }
-   
-      if (jetVectnew[index].DeltaR(jetVect[0]) > fMinDistance) return;
-    }
-    //Let's add some taste to jet and simulate pt of charged alone 
-    //eta and phi are kept as original
-    //Play a Normal Distribution
-    Float_t random = 1.;  
-    if (fSimulateChJetPt){
-      while(1){
-        random = gRandom->Gaus(0.6,0.25);
-        if (random > 0. && random < 1. && 
-            (random * jetVectnew[index].Pt()>6.)) break;
-      }
-    }
-    
-    //Set new Pt & Fill histogram accordingly
-    maxPtJet1 = random * jetVectnew[index].Pt();
-    
-    
-    fhEleadingPt->Fill( maxPtJet1 );
-
-    if (fUseAliStack){//Try Stack Information to perform UE analysis
-    
-      AliStack* mcStack = mcEvent->Stack();//Load Stack
-      Int_t nTracksMC = mcStack->GetNtrack();
-      for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
-        //Cuts
-        if(!(mcStack->IsPhysicalPrimary(iTracks))) continue;
-        
-        TParticle* mctrk = mcStack->Particle(iTracks);
-        
-        Double_t charge = mctrk->GetPDG()->Charge();
-        if (charge == 0) continue;
-        
-        if ( fUseSingleCharge ) { // Charge selection
-          if ( fUsePositiveCharge && charge < 0.) continue; // keep Positives
-          if ( !fUsePositiveCharge && charge > 0.) continue; // keep Negatives
-        }
-        
-        //Kinematics cuts on particle
-        if ((mctrk->Pt() < fTrackPtCut) || (TMath::Abs(mctrk->Eta()) > fTrackEtaCut )) continue;
-        
-        Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 ||
-                          TMath::Abs(mctrk->GetPdgCode())==2212 ||
-                          TMath::Abs(mctrk->GetPdgCode())==321;
-        
-        if ( fUseChargeHadrons && !isHadron ) continue;
-        
-        TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
+    fAnaUE->AnalyseMC(jetVect,mcEvent,pythiaGenHeader, fConePosition, fUseAliStack, fConstrainDistance, fMinDistance);
 
-        Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
-        if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
-        fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
-        
-        fhFullRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        
-        //We are not interested on stack organization but don't loose track of info
-        TVector3 tempVector =  jetVectnew[0];
-        jetVectnew[0] = jetVectnew[index];
-        jetVectnew[index] = tempVector;
-        
-        Int_t region = IsTrackInsideRegion( jetVectnew, &partVect );  
-        
-        if (region == 1) {
-          if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
-          sumPtRegionPosit += mctrk->Pt();
-          nTrackRegionPosit++;
-          fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        }
-        if (region == -1) {
-          if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
-          sumPtRegionNegat += mctrk->Pt();
-          nTrackRegionNegat++;
-          fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        }
-        if (region == 2){ //forward
-          sumPtRegionForward += mctrk->Pt();
-          nTrackRegionForward++;
-          fhRegForwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        }
-        if (region == -2){ //backward
-          sumPtRegionBackward += mctrk->Pt();
-          nTrackRegionBackward++;
-          fhRegBackwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        }
-        
-      }  // end loop stack Particles
-      
-    }else{//Try mc Particle
-
-      TClonesArray* farray = (TClonesArray*)fAOD->FindListObject("mcparticles");
-       
-      Int_t ntrks = farray->GetEntries();
-      if (fDebug>1) AliInfo(Form("In UE MC analysis tracks %d \n",ntrks));
-      for(Int_t i =0 ; i < ntrks; i++){   
-        AliAODMCParticle* mctrk = (AliAODMCParticle*)farray->At(i);
-        //Cuts
-        if (!(mctrk->IsPhysicalPrimary())) continue;
-        //if (!(mctrk->IsPrimary())) continue;
-        
-        if (mctrk->Charge() == 0 || mctrk->Charge()==-99) continue;
-        
-        if (mctrk->Pt() < fTrackPtCut ) continue;
-        if( TMath::Abs(mctrk->Eta()) > fTrackEtaCut ) continue;
-        
-        Bool_t isHadron = TMath::Abs(mctrk->GetPdgCode())==211 ||
-                          TMath::Abs(mctrk->GetPdgCode())==2212 ||
-                          TMath::Abs(mctrk->GetPdgCode())==321;
-        
-        if ( fUseChargeHadrons && !isHadron ) continue;
-        
-        TVector3 partVect(mctrk->Px(), mctrk->Py(), mctrk->Pz());
-
-        Double_t deltaPhi = jetVectnew[index].DeltaPhi(partVect)+k270rad;
-        if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
-        fhdNdEtaPhiDist->Fill( deltaPhi, maxPtJet1 );
-
-        fhFullRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        
-        //We are not interested on stack organization but don't loose track of info
-        TVector3 tempVector =  jetVectnew[0];
-        jetVectnew[0] = jetVectnew[index];
-        jetVectnew[index] = tempVector;
-        
-        Int_t region = IsTrackInsideRegion( jetVectnew, &partVect );  
-        
-        if (region == 1) { //right
-          if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
-          sumPtRegionPosit += mctrk->Pt();
-          nTrackRegionPosit++;
-          fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        }
-        if (region == -1) { //left
-          if( maxPartPtRegion < mctrk->Pt() ) maxPartPtRegion = mctrk->Pt();
-          sumPtRegionNegat += mctrk->Pt();
-          nTrackRegionNegat++;
-          fhTransRegPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        }
-        if (region == 2){ //forward
-          sumPtRegionForward += mctrk->Pt();
-          nTrackRegionForward++;
-          fhRegForwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        }
-        if (region == -2){ //backward
-          sumPtRegionBackward += mctrk->Pt();
-          nTrackRegionBackward++;
-          fhRegBackwardPartPtDistVsEt->Fill( mctrk->Pt(), maxPtJet1 );
-        }
-        
-      }//end loop AliAODMCParticle tracks
-    }
-  }
-  
-  Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.;
-  Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.;
-  if( avePosRegion > aveNegRegion ) {
-     FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg );
-  } else {
-     FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg );
-  }
-  
-  //How quantities will be sorted before Fill Min and Max Histogram
-  //  1=Plots will be CDF-like
-  //  2=Plots will be Marchesini-like
-  //  3=Minimum zone is selected as the one having lowest pt per track 
-  if( fOrdering == 1 ) {
-    if( sumPtRegionPosit > sumPtRegionNegat ) {
-      FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
-    } else {
-      FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
-    }
-    if (nTrackRegionPosit > nTrackRegionNegat ) {
-      FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
-    } else {
-      FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
-    }
-  } else if( fOrdering == 2 ) {
-    if (sumPtRegionPosit > sumPtRegionNegat) {
-      FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
-      FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
-    } else {
-      FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
-      FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
-    }
-  } else if( fOrdering == 3 ){
-     if (avePosRegion > aveNegRegion) {
-        FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg );
-        FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg );
-     }else{
-        FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg );
-        FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg );
-     }
   }
-  fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
-  
-  // Compute pedestal like magnitudes
-  fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/(2.0*fAreaReg));
-  fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/(2.0*fAreaReg));
-  // Transverse as a whole
-  fhRegTransMult->Fill( maxPtJet1, nTrackRegionPosit + nTrackRegionNegat, (nTrackRegionPosit + nTrackRegionNegat)/(2.0*fAreaReg));
-  fhRegTransSumPtVsMult->Fill(maxPtJet1, nTrackRegionPosit + nTrackRegionNegat , (sumPtRegionNegat + sumPtRegionPosit)/(2.0*fAreaReg) );
-  
-  // Fill Histograms for Forward and away side w.r.t. leading jet direction
-  // Pt dependence
-  //fhRegForwardSumPtVsEt->Fill( maxPtJet1, sumPtRegionForward/normArea );
-  //fhRegForwardMultVsEt->Fill( maxPtJet1, nTrackRegionForward/normArea );
-  //fhRegBackwardSumPtVsEt->Fill( maxPtJet1, sumPtRegionBackward/normArea );
-  //fhRegBackwardMultVsEt->Fill( maxPtJet1, nTrackRegionBackward/normArea );
-  // Multiplicity dependence
-  fhRegForwardMult->Fill(maxPtJet1, nTrackRegionForward, nTrackRegionForward/normArea);
-  fhRegForwardSumPtvsMult->Fill(maxPtJet1, nTrackRegionForward,sumPtRegionForward/normArea);
-  fhRegBackwardMult->Fill(maxPtJet1, nTrackRegionBackward, nTrackRegionBackward/normArea );
-  fhRegBackwardSumPtvsMult->Fill(maxPtJet1, nTrackRegionBackward,sumPtRegionBackward/normArea);
-  
-}
-
-//____________________________________________________________________
-void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
-{
-  // Fill sumPt of control regions
-  
-  // Max cone
-  fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
-  // Min cone
-  fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
-  // MAke distributions for UE comparison with MB data
-  fhMinRegSumPt->Fill(ptMin);
-  fhMinRegSumPtJetPtBin->Fill(leadingE, ptMin);
-  fhMaxRegSumPtJetPtBin->Fill(leadingE, ptMax);
-}
 
-//____________________________________________________________________
-void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
-{
-  // Fill average particle Pt of control regions
-  
-  // Max cone
-  fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax );
-  // Min cone
-  fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin );
-  // MAke distributions for UE comparison with MB data
-  fhMinRegAvePt->Fill(ptMin);
-}
-
-//____________________________________________________________________
-void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
-{
-  // Fill Nch multiplicity of control regions
+  fAnaUE->FillRegions(fIsNorm2Area, jetVect);
   
-  // Max cone
-  fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
-  fhRegionMultMax->Fill( nTrackPtmax );
-  // Min cone
-  fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
-  fhRegionMultMin->Fill( nTrackPtmin );
-  // MAke distributions for UE comparison with MB data
-  fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin);
 }
 
 //____________________________________________________________________
-Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect) 
-{  
-  // return de region in delta phi
-  // -1 negative delta phi 
-  //  1 positive delta phi
-  //  0 outside region
-  static const Double_t k60rad  = 60.*TMath::Pi()/180.;
-  static const Double_t k120rad = 120.*TMath::Pi()/180.;
-  
-  Int_t region = 0;
-  if( fRegionType == 1 ) {
-    if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
-    // transverse regions
-    if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; //left
-    if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;    //right
-    if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) < k60rad ) region = 2;    //forward
-    if (TMath::Abs(jetVect[0].DeltaPhi(*partVect)) > k120rad ) region = -2; //backward
-    
-  } else if( fRegionType == 2 ) {
-    // Cone regions
-    Double_t deltaR = 0.;
-    
-    TVector3 positVect,negatVect;
-    if (fConePosition==1){
-      positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
-      negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
-    }else if (fConePosition==2){
-       if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
-       positVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()+TMath::PiOver2());
-       negatVect.SetMagThetaPhi(1, 2.*atan(exp(-(jetVect[0].Eta()+jetVect[1].Eta())/2.)), jetVect[0].Phi()-TMath::PiOver2());
-    }else if (fConePosition==3){
-       if(fAnaType<2) AliFatal("Prevent error in Analysis type there might be only 1 jet. To avoid overflow better Correct UE config");
-       Double_t weightEta = jetVect[0].Eta() * jetVect[0].Pt()/(jetVect[0].Pt() + jetVect[1].Pt()) + 
-                            jetVect[1].Eta() * jetVect[1].Pt()/(jetVect[0].Pt() + jetVect[1].Pt());
-       //Double_t weightEta = jetVect[0].Eta() * jetVect[0].Mag()/(jetVect[0].Mag() + jetVect[1].Mag()) + 
-       //                     jetVect[1].Eta() * jetVect[1].Mag()/(jetVect[0].Mag() + jetVect[1].Mag());
-       positVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()+TMath::PiOver2());
-       negatVect.SetMagThetaPhi(1, 2.*atan(exp(-weightEta)), jetVect[0].Phi()-TMath::PiOver2());
-    }
-    if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) { 
-      region = 1;  
-      deltaR = positVect.DrEtaPhi(*partVect);
-    } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) { 
-      region = -1;  
-      deltaR = negatVect.DrEtaPhi(*partVect);
-    }
-    
-    if (deltaR > fConeRadius) region = 0;
-    
-  } else 
-    AliError("Unknow region type");
-  
-  // For debug (to be removed)
-  if( fDebug > 5 && region != 0 ) AliInfo(Form("Delta Phi = %3.2f region = %d \n", jetVect[0].DeltaPhi(*partVect), region)); //fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) );
+AliAnalysisTaskUE* AliAnalysisTaskUE::Instance()
+{ 
   
-  return region;
-}
-
-//____________________________________________________________________
-TObjArray*  AliAnalysisTaskUE::SortChargedParticles()
-{
-  //  return an array with all charged particles ordered according to their pT .
-  Int_t nTracks = fAOD->GetNTracks();
-  if( !nTracks ) return 0;
-  TObjArray* tracks = new TObjArray(nTracks);
-
-  for (Int_t ipart=0; ipart<nTracks; ++ipart) {
-    AliAODTrack* part = fAOD->GetTrack( ipart );
-    if( !part->TestFilterBit( fFilterBit ) ) continue; // track cut selection
-    if( !part->Charge() ) continue;
-    if( part->Pt() < fTrackPtCut ) continue;
-    tracks->AddLast( part );
-  }
-  QSortTracks( *tracks, 0, tracks->GetEntriesFast() );
-
-  nTracks = tracks->GetEntriesFast();
-  if( !nTracks ) return 0;
-
-  return tracks;
+  //Create instance
+  if (fgTaskUE) {
+       return fgTaskUE;
+  } else {
+       fgTaskUE = new AliAnalysisTaskUE();
+       return fgTaskUE;
+       }
 }
 
 //____________________________________________________________________
-TObjArray*  AliAnalysisTaskUE::FindChargedParticleJets()
+void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
 {
-  // Return a TObjArray of "charged particle jets"
-  //
-  // Charged particle jet definition from reference:
-  //
-  // "Charged jet evolution and the underlying event
-  //  in proton-antiproton collisions at 1.8 TeV"
-  //  PHYSICAL REVIEW D 65 092002, CDF Collaboration
-  //
-  // We defined "jets" as circular regions in eta-phi space with
-  // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ).
-  // Our jet algorithm is as follows:
-  //   1- Order all charged particles according to their pT .
-  //   2- Start with the highest pT particle and include in the jet all
-  //      particles within the radius R=0.7 considering each particle
-  //      in the order of decreasing pT and recalculating the centroid
-  //      of the jet after each new particle is added to the jet .
-  //   3- Go to the next highest pT particle not already included in
-  //      a jet and add to the jet all particles not already included in
-  //      a jet within R=0.7.
-  //   4- Continue until all particles are in a jet.
-  // We defined the transverse momentum of the jet to be
-  // the scalar pT sum of all the particles within the jet, where pT
-  // is measured with respect to the beam axis
-  
-  //  1 - Order all charged particles according to their pT .
-  Int_t nTracks = fAOD->GetNTracks();
-  if( !nTracks ) return 0;
-  TObjArray tracks(nTracks);
-  
-  for (Int_t ipart=0; ipart<nTracks; ++ipart) {
-    AliAODTrack* part = fAOD->GetTrack( ipart );
-    if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection
-    if( !part->Charge() ) continue;
-    if( part->Pt() < fTrackPtCut ) continue;
-    tracks.AddLast(part);
-  }
-  QSortTracks( tracks, 0, tracks.GetEntriesFast() );
-  
-  nTracks = tracks.GetEntriesFast();
-  if( !nTracks ) return 0;
-
-  TObjArray *jets = new TObjArray(nTracks);
-  TIter itrack(&tracks);
-  while( nTracks ) {
-    // 2- Start with the highest pT particle ...
-    Float_t px,py,pz,pt; 
-    AliAODTrack* track = (AliAODTrack*)itrack.Next();
-    if( !track ) continue;
-    px = track->Px();
-    py = track->Py();
-    pz = track->Pz();
-    pt = track->Pt(); // Use the energy member to store Pt
-    jets->AddLast( new TLorentzVector(px, py, pz, pt) );
-    tracks.Remove( track );
-    TLorentzVector* jet = (TLorentzVector*)jets->Last();
-    jet->SetPtEtaPhiE( 1., jet->Eta(), jet->Phi(), pt );
-    // 3- Go to the next highest pT particle not already included...
-    AliAODTrack* track1;
-    Double_t fPt = jet->E();
-    while ( (track1  = (AliAODTrack*)(itrack.Next())) ) {
-      Double_t tphi = track1->Phi(); // here Phi is from 0 <-> 2Pi
-      if (tphi > TMath::Pi()) tphi -= 2. * TMath::Pi(); // convert to  -Pi <-> Pi
-      Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-tphi);
-      Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) +
-                               dphi*dphi );
-      if( r < fConeRadius ) {
-        fPt   = jet->E()+track1->Pt();  // Scalar sum of Pt
-        // recalculating the centroid
-        Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()*track1->Pt()/fPt;
-        Double_t phi = jet->Phi()*jet->E()/fPt + tphi*track1->Pt()/fPt;
-        jet->SetPtEtaPhiE( 1., eta, phi, fPt );
-        tracks.Remove( track1 );
-      }
-    }
-    
-    tracks.Compress();
-    nTracks = tracks.GetEntries();
-    //   4- Continue until all particles are in a jet.
-    itrack.Reset();
-  } // end while nTracks
-  
-  // Convert to AODjets....
-  Int_t njets = jets->GetEntriesFast();
-  TObjArray* aodjets = new TObjArray(njets);
-  aodjets->SetOwner(kTRUE);
-  for(Int_t ijet=0; ijet<njets; ++ijet) {
-    TLorentzVector* jet = (TLorentzVector*)jets->At(ijet);
-    if (jet->E() < fChJetPtMin) continue;
-    Float_t px, py,pz,en; // convert to 4-vector
-    px = jet->E() * TMath::Cos(jet->Phi());  // Pt * cos(phi)
-    py = jet->E() * TMath::Sin(jet->Phi());  // Pt * sin(phi)
-    pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta())));
-    en = TMath::Sqrt(px * px + py * py + pz * pz);
-
-    aodjets->AddLast( new AliAODJet(px, py, pz, en) );
-  }
-  jets->Delete();
-  delete jets;
   
-  // Order jets according to their pT .
-  QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() );
-  
-  // debug
-  if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets));
-  
-  return aodjets;
-}
+  // Terminate analysis
 
-//____________________________________________________________________
-void  AliAnalysisTaskUE::QSortTracks(TObjArray &a, Int_t first, Int_t last)
-{
-  // Sort array of TObjArray of tracks by Pt using a quicksort algorithm.
-  
-  static TObject *tmp;
-  static int i;           // "static" to save stack space
-  int j;
-  
-  while (last - first > 1) {
-    i = first;
-    j = last;
-    for (;;) {
-      while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() )
-        ;
-      while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() )
-        ;
-      if (i >= j)
-        break;
-      
-      tmp  = a[i];
-      a[i] = a[j];
-      a[j] = tmp;
-    }
-    if (j == first) {
-      ++first;
-      continue;
-    }
-    tmp = a[first];
-    a[first] = a[j];
-    a[j] = tmp;
-    if (j - first < last - (j + 1)) {
-      QSortTracks(a, first, j);
-      first = j + 1;   // QSortTracks(j + 1, last);
-    } else {
-      QSortTracks(a, j + 1, last);
-      last = j;        // QSortTracks(first, j);
-    }
+  fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
+  if (!fListOfHistos){
+       AliError("Histogram List is not available");
+       return;
   }
-}
-
-//____________________________________________________________________
-void AliAnalysisTaskUE::SetRegionArea(TVector3 *jetVect)
-{
-  Double_t fAreaCorrFactor=0.;
-  Double_t deltaEta = 0.;
-  if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.;
-  else if (fRegionType==2){ 
-    deltaEta = 0.9-TMath::Abs(jetVect[0].Eta());
-    if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius;
-    else{
-      fAreaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
-      (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
-      fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-fAreaCorrFactor;
-    }
-  }else AliWarning("Unknown Rgion Type");
-  if (fDebug>10) AliInfo(Form("\n dEta=%5.3f Angle =%5.3f Region Area = %5.3f Corr Factor=%5.4f \n",deltaEta,TMath::ACos(deltaEta/fConeRadius),fAreaReg,fAreaCorrFactor));
-}
-
-//____________________________________________________________________
-void  AliAnalysisTaskUE::CreateHistos()
-{
-  fListOfHistos = new TList();
-  
-  
-  fhNJets = new TH1F("fhNJets", "n Jet",  20, 0, 20);
-  fhNJets->SetXTitle("# of jets");
-  fhNJets->Sumw2();
-  fListOfHistos->Add( fhNJets );                 // At(0) 
-  
-  fhEleadingPt  = new TH1F("hEleadingPt",   "Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhEleadingPt->SetXTitle("P_{T} (GeV/c)");
-  fhEleadingPt->SetYTitle("dN/dP_{T}");
-  fhEleadingPt->Sumw2();
-  fListOfHistos->Add( fhEleadingPt );            // At(1)
-  
-  fhMinRegPtDist = new TH1F("hMinRegPtDist",   "P_{T} distribution in Min zone",  50,0.,20.);
-  fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)");
-  fhMinRegPtDist->SetYTitle("dN/dP_{T}");
-  fhMinRegPtDist->Sumw2();
-  fListOfHistos->Add( fhMinRegPtDist );          // At(2)
-  
-  fhRegionMultMin = new TH1F("hRegionMultMin",      "N_{ch}^{90, min}",  21, -0.5,   20.5);
-  fhRegionMultMin->SetXTitle("N_{ch tracks}");
-  fhRegionMultMin->Sumw2();
-  fListOfHistos->Add( fhRegionMultMin );         // At(3)            
-  
-  fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT",  50, 0.,   20.);
-  fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)");
-  fhMinRegAvePt->Sumw2();
-  fListOfHistos->Add( fhMinRegAvePt );           // At(4)
-  
-  fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ",  50, 0.,   20.);
-  fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
-  fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
-  fhMinRegSumPt->Sumw2();
-  fListOfHistos->Add( fhMinRegSumPt );           // At(5)
-  
-  fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ",  50, 0.,   20.);
-  fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
-  fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
-  fhMinRegMaxPtPart->Sumw2();
-  fListOfHistos->Add( fhMinRegMaxPtPart );       // At(6)
-  
-  fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ",  21, -0.5,   20.5);
-  fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
-  fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
-  fhMinRegSumPtvsMult->Sumw2();
-  fListOfHistos->Add( fhMinRegSumPtvsMult );     // At(7);
-  
-  fhdNdEtaPhiDist  = new TH2F("hdNdEtaPhiDist", Form("Charge particle density |#eta|<%3.1f vs #Delta#phi", fTrackEtaCut),  
-                               62, 0.,   2.*TMath::Pi(), fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
-  fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
-  fhdNdEtaPhiDist->SetYTitle("Leading Jet P_{T}");
-  fhdNdEtaPhiDist->Sumw2();
-  fListOfHistos->Add( fhdNdEtaPhiDist );        // At(8)
-  
-  // Can be use to get part pt distribution for differente Jet Pt bins
-  fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
-                                     200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
-  fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
-  fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
-  fhFullRegPartPtDistVsEt->Sumw2();
-  fListOfHistos->Add( fhFullRegPartPtDistVsEt );  // At(9) 
-  
-  // Can be use to get part pt distribution for differente Jet Pt bins
-  fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", Form( "dN/dP_{T} in tranvese regions |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),  
-                                     200,0.,100., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
-  fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
-  fhTransRegPartPtDistVsEt->Sumw2();
-  fListOfHistos->Add( fhTransRegPartPtDistVsEt );  // At(10)
-  
-  
-  fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt",  "P_{T}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegionSumPtMaxVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionSumPtMaxVsEt );      // At(11)
-  
-  fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt",   "P_{T}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegionSumPtMinVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionSumPtMinVsEt );      // At(12)
-  
-  fhRegionMultMax = new TH1I("hRegionMultMax",      "N_{ch}^{90, max}",  21, -0.5,   20.5);
-  fhRegionMultMax->SetXTitle("N_{ch tracks}");
-  fhRegionMultMax->Sumw2();
-  fListOfHistos->Add( fhRegionMultMax );           // At(13)
-  
-  fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt",  "N_{ch}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
-  fhRegionMultMaxVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionMultMaxVsEt );       // At(14)
-  
-  fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt",  "N_{ch}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
-  fhRegionMultMinVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionMultMinVsEt );      // At(15)
-  
-  fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegionAveSumPtVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionAveSumPtVsEt );     // At(16)
-  
-  fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegionDiffSumPtVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionDiffSumPtVsEt );    // At(17)
-  
-  fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegionAvePartPtMaxVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionAvePartPtMaxVsEt );  // At(18)
-  
-  fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegionAvePartPtMinVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionAvePartPtMinVsEt );   // At(19)
-  
-  fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegionMaxPartPtMaxVsEt->Sumw2();
-  fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt );    // At(20)
-  /*
-  fhRegForwardSumPtVsEt = new TH1F("hRegForwardSumPtVsEt", "Forward #sum{p_{T}} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegForwardSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegForwardSumPtVsEt->Sumw2();
-  fListOfHistos->Add( fhRegForwardSumPtVsEt );    // At(21)
-  
-  fhRegForwardMultVsEt = new TH1F("hRegForwardMultVsEt", "Forward #sum{N_{ch}} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegForwardMultVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegForwardMultVsEt->Sumw2();
-  fListOfHistos->Add( fhRegForwardMultVsEt );    // At(22)
-  
-  fhRegBackwardSumPtVsEt = new TH1F("hRegBackwardSumPtVsEt", "Backward #sum{p_{T}} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegBackwardSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegBackwardSumPtVsEt->Sumw2();
-  fListOfHistos->Add( fhRegBackwardSumPtVsEt );    // At(23)
-  
-  fhRegBackwardMultVsEt = new TH1F("hRegBackwardMultVsEt", "Backward #sum{N_{ch}} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-  fhRegBackwardMultVsEt->SetXTitle("P_{T} (GeV/c)");
-  fhRegBackwardMultVsEt->Sumw2();
-  fListOfHistos->Add( fhRegBackwardMultVsEt );    // At(24)
-                                  */
-  fhRegForwardMult = new TH2F("hRegForwardMult",      "N_{ch}^{forward}",  
-                              fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
-  fhRegForwardMult->SetXTitle("N_{ch tracks}");
-  fhRegForwardMult->Sumw2();
-  fListOfHistos->Add( fhRegForwardMult );           // At(25)
-  
-  fhRegForwardSumPtvsMult = new TH2F("hRegForwardSumPtvsMult", "Forward #Sigma p_{T} vs. Multiplicity ",  
-                                     fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
-  fhRegForwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
-  fhRegForwardSumPtvsMult->SetXTitle("N_{charge}");
-  fhRegForwardSumPtvsMult->Sumw2();
-  fListOfHistos->Add( fhRegForwardSumPtvsMult );     // At(26);
-  
-  fhRegBackwardMult = new TH2F("hRegBackwardMult",      "N_{ch}^{backward}",  
-                               fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
-  fhRegBackwardMult->SetXTitle("N_{ch tracks}");
-  fhRegBackwardMult->Sumw2();
-  fListOfHistos->Add( fhRegBackwardMult );           // At(27)
-  
-  fhRegBackwardSumPtvsMult = new TH2F("hRegBackwardSumPtvsMult", "Backward #Sigma p_{T} vs. Multiplicity ",  
-                                      fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
-  fhRegBackwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
-  fhRegBackwardSumPtvsMult->SetXTitle("N_{charge}");
-  fhRegBackwardSumPtvsMult->Sumw2();
-  fListOfHistos->Add( fhRegBackwardSumPtvsMult );     // At(28);
-  
-  fhRegForwardPartPtDistVsEt = new TH2F("hRegForwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
-                                       200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
-  fhRegForwardPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
-  fhRegForwardPartPtDistVsEt->SetXTitle("p_{T}");
-  fhRegForwardPartPtDistVsEt->Sumw2();
-  fListOfHistos->Add( fhRegForwardPartPtDistVsEt );  // At(29) 
-  
-  fhRegBackwardPartPtDistVsEt = new TH2F("hRegBackwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading Jet P_{T}", fTrackEtaCut),
-                                         200,0.,100., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist);
-  fhRegBackwardPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
-  fhRegBackwardPartPtDistVsEt->SetXTitle("p_{T}");
-  fhRegBackwardPartPtDistVsEt->Sumw2();
-  fListOfHistos->Add( fhRegBackwardPartPtDistVsEt );  // At(30) 
-  
-  fhRegTransMult  = new TH2F("hRegTransMult",      "N_{ch}^{transv}",  
-                             fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
-  fhRegTransMult->SetXTitle("N_{ch tracks}");
-  fhRegTransMult->Sumw2();
-  fListOfHistos->Add( fhRegTransMult );           // At(31)
-  
-  fhRegTransSumPtVsMult = new TH2F("hRegTransSumPtVsMult", "Transverse #Sigma p_{T} vs. Multiplicity ",  
-                                   fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 21, -0.5,   20.5);
-  fhRegTransSumPtVsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
-  fhRegTransSumPtVsMult->SetXTitle("N_{charge}");
-  fhRegTransSumPtVsMult->Sumw2();
-  fListOfHistos->Add( fhRegTransSumPtVsMult );     // At(32);
-  
-  fhMinRegSumPtJetPtBin = new TH2F("hMinRegSumPtJetPtBin",      "Transverse Min Reg #Sigma p_{T} per jet Pt Bin",  
-                                   fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 50, 0.,   20.);
-  fhMinRegSumPtJetPtBin->SetXTitle("Leading Jet P_{T}");
-  fhMinRegSumPtJetPtBin->Sumw2();
-  fListOfHistos->Add( fhMinRegSumPtJetPtBin );           // At(33)
-  
-  fhMaxRegSumPtJetPtBin = new TH2F("hMaxRegSumPtJetPtBin",      "Transverse Max Reg #Sigma p_{T} per jet Pt Bin",  
-                                   fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, 50, 0.,   20.);
-  fhMaxRegSumPtJetPtBin->SetXTitle("Leading Jet P_{T}");
-  fhMaxRegSumPtJetPtBin->Sumw2();
-  fListOfHistos->Add( fhMaxRegSumPtJetPtBin );           // At(34)
-  
-  fhVertexMult = new TH1F("hVertexMult",      "Multiplicity in Main Vertex", 81, -0.5 , 80.5);
-  fhVertexMult->SetXTitle("Main Vertex Multiplicity");
-  fhVertexMult->Sumw2();
-  fListOfHistos->Add( fhVertexMult ); //At(35)
-  
-  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); 
-  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
-  fh1Xsec->Sumw2();
-  fListOfHistos->Add( fh1Xsec );            //At(36)
-  
-  fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
-  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
-  fh1Trials->Sumw2();
-  fListOfHistos->Add( fh1Trials ); //At(37)
-  
-  fSettingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
-  fSettingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
-  fSettingsTree->Branch("fTrigger", &fTrigger,"TriggerFlag/I");
-  fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
-  fSettingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
-  fSettingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
-  fSettingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
-  fSettingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
-  fSettingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
-  fSettingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
-  fSettingsTree->Branch("fAnaType", &fAnaType, "Ana/I");        
-  fSettingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
-  fSettingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
-  fSettingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
-  fSettingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
-  fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
-  fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
-  fSettingsTree->Fill();
-
-  
-  fListOfHistos->Add( fSettingsTree );    // At(38)
-  
-  /*   
-   // For debug region selection
-   fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi",      
-   20, -2.,2., 62, -TMath::Pi(),   TMath::Pi());
-   fhValidRegion->SetYTitle("#Delta#phi");
-   fhValidRegion->Sumw2();
-   fListOfHistos->Add( fhValidRegion );  // At(15)
-   */
-}
-
-
-
-//____________________________________________________________________
-void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
-{
-  // Terminate analysis
-  //
-  
-  //Save Analysis Settings
-  //  WriteSettings();
-  
-  // Normalize histos to region area TODO: 
-  // Normalization done at Analysis, taking into account 
-  // area variations on per-event basis (cone case)
-   
-   //HIGH WARNING!!!!!: DO NOT SCALE ANY OF THE ORIGINAL HISTOGRAMS
-   //MAKE A COPY, DRAW IT, And later sacale that copy. CAF Issue!!!!!
-  
-  // Draw some Test plot to the screen
-  if (!gROOT->IsBatch()) {
-    
-    // Update pointers reading them from the output slot
-    fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
-    if( !fListOfHistos ) {
-      AliError("Histogram List is not available");
-      return;
-    }
-    fhNJets              = (TH1F*)fListOfHistos->At(0);
-    fhEleadingPt         = (TH1F*)fListOfHistos->At(1);
-    fhdNdEtaPhiDist      = (TH2F*)fListOfHistos->At(8);
-    fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11);
-    fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12);
-    fhRegionMultMaxVsEt  = (TH1F*)fListOfHistos->At(14);
-    fhRegionMultMinVsEt  = (TH1F*)fListOfHistos->At(15);
-    fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16);
-    //fhRegForwardSumPtVsEt = (TH1F*)fListOfHistos->At(21);
-    //fhRegBackwardSumPtVsEt = (TH1F*)fListOfHistos->At(23);
-    
-    //fhValidRegion  = (TH2F*)fListOfHistos->At(21);
-    
-    // Canvas
-    TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1100,700);
-    c1->Divide(2,2);
-    c1->cd(1);
-    TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
-    h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
-    //h1r->Scale( areafactor );
-    h1r->SetMarkerStyle(20);
-    h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
-    h1r->SetYTitle("P_{T}^{90, max}");
-    h1r->DrawCopy("p");
-    
-    c1->cd(2);
-    
-    TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
-    h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
-    //h2r->Scale( areafactor );
-    h2r->SetMarkerStyle(20);
-    h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
-    h2r->SetYTitle("P_{T}^{90, min}");
-    h2r->DrawCopy("p");
-    
-    c1->cd(3);
-    TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
-    //TH1F *h41r = new TH1F("hRegForwvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
-    //TH1F *h42r = new TH1F("hRegBackvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
-    //h41r->Divide(fhRegForwardSumPtVsEt,fhEleadingPt,1,1);
-    //h42r->Divide(fhRegBackwardSumPtVsEt,fhEleadingPt,1,1);
-    h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
-    //h4r->Scale(2.); // make average
-    //h4r->Scale( areafactor );
-    h4r->SetYTitle("#DeltaP_{T}^{90}");
-    h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
-    h4r->SetMarkerStyle(20);
-    h4r->DrawCopy("p");
-    /*h41r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
-    h41r->SetMarkerStyle(22);
-    h41r->DrawCopy("p same");
-    h42r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
-    h42r->SetMarkerStyle(23);
-    h42r->SetMarkerColor(kRed);
-    h42r->DrawCopy("p same");
-    */
-    c1->cd(4);
-    TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-    TH1F *h6r = new TH1F("hRegionMultMinVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
-    
-    h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
-    h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
-    //h5r->Scale( areafactor );
-    //h6r->Scale( areafactor );
-    h5r->SetYTitle("N_{Tracks}^{90}");
-    h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
-    h5r->SetMarkerStyle(20);
-    h5r->DrawCopy("p");
-    h6r->SetMarkerStyle(21);
-    h6r->SetMarkerColor(2);
-    h6r->SetYTitle("N_{Tracks}^{90}");
-    h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
-    h6r->DrawCopy("p same");
-    c1->Update();
-    
-    //Get Normalization
-    fh1Xsec                = (TProfile*)fListOfHistos->At(21);
-    fh1Trials              = (TH1F*)fListOfHistos->At(22);
-    
-    Double_t xsec = fh1Xsec->GetBinContent(1);
-    Double_t ntrials = fh1Trials->GetBinContent(1);
-    Double_t normFactor = xsec/ntrials;
-    if(fDebug > 1)Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor);
-    
-    
-    
-    TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
-    TH1 * copy = 0;
-    c2->Divide(2,2);
-    c2->cd(1);
-    fhEleadingPt->SetMarkerStyle(20);
-    fhEleadingPt->SetMarkerColor(2);
-    //if( normFactor > 0.) fhEleadingPt->Scale(normFactor);
-    //fhEleadingPt->Draw("p");
-    copy = fhEleadingPt->DrawCopy("p");
-    if( normFactor > 0.) copy->Scale(normFactor);
-    gPad->SetLogy();
-    
-    c2->cd(2);
-    Int_t xbin1 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(fMinJetPtInHist);
-    Int_t xbin2 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(fMaxJetPtInHist);
-    TH1D* dNdEtaPhiDistAllJets = fhdNdEtaPhiDist->ProjectionX("dNdEtaPhiDistAllJets",xbin1,xbin2);
-    dNdEtaPhiDistAllJets->SetMarkerStyle(20);
-    dNdEtaPhiDistAllJets->SetMarkerColor(2);
-    dNdEtaPhiDistAllJets->DrawCopy("p");
-    gPad->SetLogy();
-    
-    c2->cd(3);      
-    fhNJets->DrawCopy();
-    
-    //c2->cd(4);      
-    //fhValidRegion->DrawCopy("p");
-    
-    //fhTransRegPartPtDist  = (TH1F*)fListOfHistos->At(2);
-    fhRegionMultMin          = (TH1F*)fListOfHistos->At(3);
-    fhMinRegAvePt     = (TH1F*)fListOfHistos->At(4);
-    fhMinRegSumPt     = (TH1F*)fListOfHistos->At(5);   
-    //fhMinRegMaxPtPart     = (TH1F*)fListOfHistos->At(6);
-    fhMinRegSumPtvsMult   = (TH1F*)fListOfHistos->At(7);
-    
-    
-    // Canvas
-    TCanvas* c3 = new TCanvas("c3"," pT dist",160,160,1200,800);
-    c3->Divide(2,2);
-    c3->cd(1);
-    /*fhTransRegPartPtDist->SetMarkerStyle(20);
-     fhTransRegPartPtDist->SetMarkerColor(2); 
-     fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
-     fhTransRegPartPtDist->DrawCopy("p");
-     gPad->SetLogy();
-     */
-    c3->cd(2); 
-    fhMinRegSumPt->SetMarkerStyle(20);
-    fhMinRegSumPt->SetMarkerColor(2);  
-    //fhMinRegSumPt->Scale(areafactor);
-    fhMinRegSumPt->DrawCopy("p");
-    gPad->SetLogy();
-    
-    c3->cd(3);
-    fhMinRegAvePt->SetMarkerStyle(20);
-    fhMinRegAvePt->SetMarkerColor(2);  
-    //fhMinRegAvePt->Scale(areafactor);
-    fhMinRegAvePt->DrawCopy("p");
-    gPad->SetLogy();
-    
-    c3->cd(4);
-    
-    TH1F *h7r = new TH1F("hRegionMultMinVsMult",   "",  21, -0.5,   20.5);
-    
-    h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
-    
-    h7r->SetMarkerStyle(20);
-    h7r->SetMarkerColor(2);   
-    h7r->DrawCopy("p");
-    
-    c3->Update();
-    
-    
-    /*    c2->cd(3);      
-     fhValidRegion->SetMarkerStyle(20);
-     fhValidRegion->SetMarkerColor(2);
-     fhValidRegion->DrawCopy("p");
-     */    
-    c2->Update();
+  if (!gROOT->IsBatch()){
+  //call class AliHistogramsUE
+  AliHistogramsUE *histos=new AliHistogramsUE(fListOfHistos);
+  histos->DrawUE(0);
   } else {
-    AliInfo(" Batch mode, not histograms will be shown...");
+        AliInfo(" Batch mode, not histograms will be shown...");
   }
-  
+
   if( fDebug > 1 ) 
     AliInfo("End analysis");
-  
 }
 
 void  AliAnalysisTaskUE::WriteSettings()
 { 
-  if (fDebug>5){
+
+  //Print analysis settings on screen
+  if (fDebug > 5){
     AliInfo(" All Analysis Settings in Saved Tree");
-    fSettingsTree->Scan();
+    fAnaUE->WriteSettings();
   }
 }
index b9a6b2c..24cd64e 100644 (file)
@@ -4,23 +4,55 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
+////////////////////////////////////////////////////////////////////////
+//
+// Analysis class for Underlying Event studies
+//
+// Look for correlations on the tranverse regions to 
+// the leading charged jet
+//
+// This class needs as input AOD with track and Jets.
+// The output is a list of histograms
+//
+// AOD can be either connected to the InputEventHandler  
+// for a chain of AOD files 
+// or 
+// to the OutputEventHandler
+// for a chain of ESD files, so this case class should be 
+// in the train after the Jet finder
+//
+//    Arian.Abrahantes.Quintana@cern.ch 
+//    Ernesto.Lopez.Torres@cern.ch
+// 
+////////////////////////////////////////////////////////////////////////
+
 #include "AliAnalysisTask.h"
 
-class AliESDEvent;
+class AliAnalyseUE;
 class AliAODEvent;
+//class AliAODHandler;
+class AliAODInputHandler;
+class AliESDEvent;
+//class AliLog;
 class TH1F;
 class TH2F;
 class TH1I;
 class TProfile;
-class TVector3;
+//class TROOT;
 class TTree;
+class TVector3;
 
 class  AliAnalysisTaskUE : public AliAnalysisTask
   {
   public:
     AliAnalysisTaskUE(const char* name="AliAnalysisTaskUE");
     virtual           ~AliAnalysisTaskUE() {if ( fListOfHistos ) delete fListOfHistos; }
-    
+    AliAnalysisTaskUE(const  AliAnalysisTaskUE &det);
+    AliAnalysisTaskUE&   operator=(const  AliAnalysisTaskUE &det);
+       
+    // return instance of the singleton
+    static  AliAnalysisTaskUE* Instance();
+      
     // Implementation of interace methods
     virtual     Bool_t Notify();
     virtual     void   ConnectInputData(Option_t *);
@@ -28,216 +60,202 @@ class  AliAnalysisTaskUE : public AliAnalysisTask
     virtual     void   Exec(Option_t *option);
     virtual     void   Terminate(Option_t *);
 
-    //Select the trigger
-    void SelectTrigger(Int_t trig) { fTrigger = trig; }
+    //  Setters/Getters
+    virtual     void         SetDebugLevel( Int_t level )  { fDebug = level; }
+    virtual const Int_t       GetDebugLevel()  { return fDebug; }
+    
+    void       SetPtRangeInHist( Int_t bin, Double_t min, Double_t max ) {
+                       fBinsPtInHist = bin; 
+                       fMinJetPtInHist = min; 
+                       fMaxJetPtInHist = max; 
+                       }
 
-    //  Setters
-    virtual     void   SetDebugLevel( Int_t level )  { fDebug = level; }
-    void   SetPtRangeInHist( Int_t bin, Double_t min, Double_t max ) {
-      fBinsPtInHist = bin; 
-      fMinJetPtInHist = min; 
-      fMaxJetPtInHist = max; 
-    }
+    // Read AODs
+    void         SelectAODBranch(const char* val)     { fAODBranch = val;   }
+    virtual const TString      GetAODBranch()                       { return fAODBranch;   }
+    virtual const AliAODEvent*   GetAOD()                             { return fAOD; } 
 
-    // Read deltaAODs
-    void   ReadDeltaAOD()                   { fDeltaAOD = kTRUE; }
-    void   SelectDeltaAODBranch(const char* val)     { fDeltaAODBranch = val;   }
-    void   SelectAODBranch(const char* val)     { fAODBranch = val;   }
+    // Setters/Getters for MC
+    void       SetConstrainDistance(Bool_t val1, Double_t val2){ fMinDistance = val2; fConstrainDistance = val1;}
+    void       SetSimulateChJetPt(){fSimulateChJetPt = kTRUE;}
+    void       SetUseAODMCParticle(){fUseAliStack = kFALSE;}
+    void       SetUseMCBranch(){fUseMCParticleBranch = kTRUE;}
+
+    virtual const Bool_t       GetConstrainDistance() {return fConstrainDistance;}
+    virtual const Double_t     GetMinDistance() {return fMinDistance;}
+    virtual const Bool_t       GetSimulateChJetPt(){return fSimulateChJetPt;}
+    virtual const Bool_t       GetUseAODMCParticle(){return fUseAliStack;}
+    virtual const Bool_t       GetUseMCParticleBranch(){return fUseMCParticleBranch;}
 
-    // Setters for MC
-    void  SetUseMCBranch(){fUseMCParticleBranch = kTRUE;}
-    void  SetConstrainDistance(Bool_t val1, Double_t val2){ fMinDistance = val2; fConstrainDistance = val1;}
-    void  SetSimulateChJetPt(){fSimulateChJetPt = kTRUE;}
-    void  SetUseAODMCParticle(){fUseAliStack = kFALSE;}
-    
     //Setters for Events QA
     void  SetZVertex( Double_t val ) { fZVertex = val; }
     void  SetTracksInVertex( Int_t val ){ fnTracksVertex = val; }
     
-    // Stters for UE Analysis
-    void   SetAnaTopology( Int_t val )    { fAnaType = val;    }
-    void   SetRegionType( Int_t val )     { fRegionType = val; }
-    void   SetUseChPartJet( Int_t val )   { fUseChPartJet = val; }
-    void   SetUseChargeHadrons( Bool_t val ){ fUseChargeHadrons = val; }
-    void   SetPtSumOrdering( Int_t val ) { fOrdering = val;   }
-    void   SetFilterBit( UInt_t val )     { fFilterBit = val;  }
-    void   SetJetsOnFly( Bool_t val )     { fJetsOnFly = val;  }
-    void   SetConeRadius( Double_t val )  { fConeRadius = val; }
-    void   SetConePosition(Int_t val)     { fConePosition= val; }
-    void   SetUseSingleCharge()  { fUseSingleCharge = kTRUE; } 
+    // Setters/Getters for UE Analysis
+    void   SetAnaTopology( Int_t val )       { fAnaType = val;    }
+    void   SetConePosition(Int_t val)        { fConePosition= val; }
+    void   SetConeRadius( Double_t val )     { fConeRadius = val; }
+    void   SetDoNotNormalizeQuantities()     { fIsNorm2Area = kFALSE; }
+    void   SetFilterBit( UInt_t val )        { fFilterBit = val;  }
+    void   SetJetsOnFly( Bool_t val )        { fJetsOnFly = val;  }
+    void   SetPtSumOrdering( Int_t val )     { fOrdering = val;   }
+    void   SetRegionType( Int_t val )        { fRegionType = val; }
+    void   SetUseChargeHadrons( Bool_t val ) { fUseChargeHadrons = val; }
+    void   SetUseChPartJet( Int_t val )      { fUseChPartJet = val; }
     void   SetUseNegativeChargeType()        { fUsePositiveCharge = kFALSE; }
-    void   SetDoNotNormalizeQuantities()  { fIsNorm2Area = kFALSE; }
+    void   SetUseSingleCharge()              { fUseSingleCharge = kTRUE; } 
+    
+    virtual const Int_t     GetAnaTopology()                    { return fAnaType; }
+    virtual const Int_t     GetConePosition()                  { return fConePosition; }
+    virtual const Double_t  GetConeRadius()                     { return fConeRadius; }
+    virtual const Bool_t    GetDoNotNormalizeQuantities()      { return fIsNorm2Area; }
+    virtual const UInt_t    GetFilterBit()                      { return fFilterBit;  }
+    virtual const Bool_t    GetJetsOnFly()                      { return fJetsOnFly;  }
+    virtual const Int_t     GetPtSumOrdering()                 { return fOrdering;   }
+    virtual const Int_t     GetRegionType()                    { return fRegionType; }
+    virtual const Bool_t    GetUseChargeHadrons()              { return fUseChargeHadrons; }
+    virtual const Int_t     GetUseChPartJet()                  { return fUseChPartJet; }
+    virtual const Bool_t    GetUseNegativeChargeType()         { return fUsePositiveCharge; }
+    virtual const Bool_t    GetUseSingleCharge()               { return fUseSingleCharge; } 
+    
     // Jet cuts
-    void   SetPtMinChPartJet( Double_t val )  { fChJetPtMin = val; }
     void   SetJet1EtaCut( Double_t val )      { fJet1EtaCut = val; }
     void   SetJet2DeltaPhiCut( Double_t val ) { fJet2DeltaPhiCut = val; }
     void   SetJet2RatioPtCut( Double_t val )  { fJet2RatioPtCut = val; }
     void   SetJet3PtCut( Double_t val )       { fJet3PtCut = val; }
+    void   SetPtMinChPartJet( Double_t val )  { fChJetPtMin = val; }
+    
+    virtual const Double_t   GetJet1EtaCut()                { return fJet1EtaCut; }
+    virtual const Double_t   GetJet2DeltaPhiCut()            { return fJet2DeltaPhiCut; }
+    virtual const Double_t   GetJet2RatioPtCut()            { return fJet2RatioPtCut; }
+    virtual const Double_t   GetJet3PtCut()                 { return fJet3PtCut; }
+    virtual const Double_t   GetPtMinChPartJet()            { return fChJetPtMin; }
+    
     // track cuts
-    void   SetTrackPtCut( Double_t val )  { fTrackPtCut = val; }
     void   SetTrackEtaCut( Double_t val ) { fTrackEtaCut = val; }
-    
+    void   SetTrackPtCut( Double_t val )  { fTrackPtCut = val; }
+
+    virtual const Double_t   GetTrackEtaCut()           { return fTrackEtaCut; }
+    virtual const Double_t   GetTrackPtCut()            { return fTrackPtCut; }
+
+  protected:
+  static AliAnalysisTaskUE*     fgTaskUE;        // Pointer to single instance
+
   private:
-    AliAnalysisTaskUE(const  AliAnalysisTaskUE &det);
-    AliAnalysisTaskUE&   operator=(const  AliAnalysisTaskUE &det);
-    void   AnalyseUE();
-    Int_t   IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect);
-    void   CreateHistos();
-    void   SetRegionArea(TVector3 *jetVect);
-    void   FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin );
-    void   FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin );
-    void   FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin );
-    TObjArray*  FindChargedParticleJets();
-    TObjArray*  SortChargedParticles();
-    void   QSortTracks(TObjArray &a, Int_t first, Int_t last);
-    void   WriteSettings();
-    
-    Int_t      fTrigger;         //Trigger flag as defined in AliAnalysisHelperJetTasks.h
-    Int_t      fDebug;           //  Debug flag
-    Bool_t      fDeltaAOD;        //  Read jets from delta AOD 
-    TString     fDeltaAODBranch;  //  Jet branch name from delta AOD
-    TString     fAODBranch;       //  Jet branch name from standard AOD
-    TClonesArray*  fArrayJets;       //!  Array of Jets from delta AOD
-
-    AliAODEvent*  fAOD;             //! AOD Event 
-    AliAODEvent*  fAODjets;         //! AOD Event for reconstructed on the fly (see ConnectInputData()
-    TList*  fListOfHistos;    //  Output list of histograms
+    
+    void         AnalyseUE();
+    void         FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin );
+    void         FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin );
+    void         FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin );
+    TObjArray*   FindChargedParticleJets();
+    Int_t        IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect);
+    void         QSortTracks(TObjArray &a, Int_t first, Int_t last);
+    void         SetRegionArea(TVector3 *jetVect);
+    TObjArray*   SortChargedParticles();
+    void         WriteSettings();
+    
+    AliAnalyseUE*  fAnaUE;           //! points to AliAnalyseUE class
+    AliAODEvent*   fAOD;             //! AOD Event 
+    TString        fAODBranch;       //  Jet branch name from standard AOD
+    Int_t          fDebug;           //  Debug flag
+
+    TList*         fListOfHistos;    //  Output list of histograms
     
     // Config
-    Int_t   fBinsPtInHist;     //  # bins for Pt histos range
-    Double_t   fMinJetPtInHist;   //  min Jet Pt value for histo range
-    Double_t   fMaxJetPtInHist;   //  max Jet Pt value for histo range
-    Bool_t     fIsNorm2Area;      // Apply Area Normalization to collected observables
+    Int_t          fBinsPtInHist;     //  # bins for Pt histos range
+    Bool_t         fIsNorm2Area;      // Apply Area Normalization to collected observables
+    Double_t       fMaxJetPtInHist;   //  max Jet Pt value for histo range
+    Double_t       fMinJetPtInHist;   //  min Jet Pt value for histo range
     
     // For MC
-    Bool_t fUseMCParticleBranch;  // Run Over mcparticles branch in AOD
-    Bool_t fConstrainDistance;    // Constrain Distance between rec jet and pyth
-    Double_t fMinDistance;  // Minimum distance between rec jet and pyth
-    Bool_t fSimulateChJetPt; // Naive simulation of charged jet Pt from original Jet in MC Header
-    Bool_t fUseAliStack;     // Use AliSatck for particle info otherwise "mcparticles" branch in AOD
+    Bool_t        fConstrainDistance;    // Constrain Distance between rec jet and pyth
+    Double_t       fMinDistance;          // Minimum distance between rec jet and pyth
+    Bool_t         fSimulateChJetPt;      // Naive simulation of charged jet Pt from original Jet in MC Header
+    Bool_t         fUseAliStack;          // Use AliSatck for particle info otherwise "mcparticles" branch in AOD
+    Bool_t         fUseMCParticleBranch;  // Run Over mcparticles branch in AOD
     
     // Cuts Events type
-    Int_t fnTracksVertex;  // QA tracks pointing to principal vertex (= 3 default) 
-    Double_t fZVertex;     // Position of Vertex in Z direction
-    
-    // Cuts 
-    Int_t   fAnaType;          // Analysis type on jet topology: 
-    //     1=inclusive  (default) 
-    //     2=back to back inclusive
-    //     3=back to back exclusive
-    //     4=Pt max (max Pt track in region)
-    //     5=gama jet (back to back) ???
-    //  Minimum bias
-    //     31 = Semi jet (charged leading particle jets)
-    //     32 = Random jetcone  ?
-    //     33 = Swiss chees   ?
+    Int_t          fnTracksVertex;        // QA tracks pointing to principal vertex (= 3 default) 
+    Double_t       fZVertex;              // Position of Vertex in Z direction
+    
+    // Cuts UE analysis
+    Int_t          fAnaType;              // Analysis type on jet topology: 
+                                               //     1=inclusive  (default) 
+                                               //     2=back to back inclusive
+                                               //     3=back to back exclusive
+                                               //     4=Pt max (max Pt track in region)
+                                               //     5=gama jet (back to back) ???
+                                               //  Minimum bias
+                                               //     31 = Semi jet (charged leading particle jets)
+                                               //     32 = Random jetcone  ?
+                                               //     33 = Swiss chees   ?
+    
+    
+    Int_t          fConePosition;         // This parameter set how will cone center in transversal zone will be set
+                                          //    1 : To be used in any jet topology (default value)
+                                          //        eta_cone = eta_leadingjet
+                                          //        phi_cone = phi_leadingjet + - 90
+                                          //    2 : To be used in multiple jet topology (code will cry otherwise)
+                                          //        eta_cone = (eta_leadingjet + eta_subleadingjet)/2
+                                          //        phi_cone = phi_leadingjet + - 90
+    
+    Double_t       fConeRadius;           // if selected Cone-like region type, set Radius (=0.7 default)
+    
+    UInt_t         fFilterBit;            // Select tracks from an specific track cut (default 0xFF all track selected)
+    
+    Bool_t         fJetsOnFly;            // if jets are reconstructed on the fly from AOD tracks (see ConnectInputData() )
     
     // UE analysis is conducted in different type of regions
     // Transverse are those like defined in: R. Field Acta Physica Polonica B. Vol 36 No. 2 pg 167 (2005) 
     // Cone regions like defined in: Phys. Rev. D 70, 072002 (2004)
-    Int_t   fRegionType;       // 1 = transverse regions (default)
-                               // 2 = cone regions   
-    Double_t   fConeRadius;       // if selected Cone-like region type, set Radius (=0.7 default)
-    Int_t   fConePosition;     // This parameter set how will cone center in transversal zone will be set
-                               //    1 : To be used in any jet topology (default value)
-                               //        eta_cone = eta_leadingjet
-                               //        phi_cone = phi_leadingjet + - 90
-                               //    2 : To be used in multiple jet topology (code will cry otherwise)
-                               //        eta_cone = (eta_leadingjet + eta_subleadingjet)/2
-                               //        phi_cone = phi_leadingjet + - 90
-    Double_t   fAreaReg;       // Area of the region To be used as normalization factor when filling histograms
-                               // if fRegionType = 2 not always it is included within eta range
-    Bool_t   fUseChPartJet;     // Use "Charged Particle Jet" instead of jets from AOD see FindChargedParticleJets()
-    Bool_t   fUseChPartMaxPt;   // Use "Charged Particle with max Pt" instead of any jets to define forward region
-
-    Bool_t   fUseChargeHadrons;   // Only use charge hadrons
+    Int_t         fRegionType;           // 1 = transverse regions (default)
+                                          // 2 = cone regions   
+    
+
+
+    Bool_t         fUseChargeHadrons;     // Only use charge hadrons
+    Bool_t         fUseChPartJet;         // Use "Charged Particle Jet" instead of jets from AOD see FindChargedParticleJets()
     
     // Theoreticians ask for tools charge-aware
     // especially those which are related to multiplicity and MC-tunings
     // see arXiv:hep-ph/0507008v3
-    Bool_t   fUseSingleCharge;     //Make analysis for a single type of charge (=kFALSE default)
-    Bool_t   fUsePositiveCharge;   //If Single type of charge used then set which one (=kTRUE default positive)
-    Int_t   fOrdering;         //  Pt and multiplicity summation ordering:
-    //     1=CDF-like -independent sorting according quantity to be scored: Double sorting- (default)
-    //       if Pt summation will be scored take Pt minimum between both zones and 
-    //          fill Pt Max. and Min. histog. accordingly
-    //       if Multiplicity summation will be scored take Mult. minimum between both zones and 
-    //          fill Mult Max and Min histog. accordingly
-    //       Bib:
-    //     2=Marchesini-like (Only Pt sorting: Single sorting)
-    //          sort only according Pt summation scored, find minimum between both zones and
-    //          fill Pt and Multiplicity Max and Min summation histog. following only this criterium
-    //       Bib: Phys. Rev. D 38, 3419 (1988)
-    //     3=Nameless pt per track single sorting
-    //          sort according to pt per track scored in each transverse zone 
-    //          lowest values indicates minimum zone.   
-    //     4=User Selection sorting (NOTE: USER must implement it within cxx)
-    
-    UInt_t   fFilterBit;        // Select tracks from an specific track cut (default 0xFF all track selected)
-    Bool_t   fJetsOnFly;        // if jets are reconstructed on the fly from AOD tracks (see ConnectInputData() )
+    Bool_t         fUsePositiveCharge;    //If Single type of charge used then set which one (=kTRUE default positive)
+    Bool_t         fUseSingleCharge;      //Make analysis for a single type of charge (=kFALSE default)
+    
+    Int_t          fOrdering;             //  Pt and multiplicity summation ordering:
+                   //     1=CDF-like -independent sorting according quantity to be scored: Double sorting- (default)
+                   //        if Pt summation will be scored take Pt minimum between both zones and 
+                   //        fill Pt Max. and Min. histog. accordingly
+                   //        if Multiplicity summation will be scored take Mult. minimum between both zones and 
+                   //        fill Mult Max and Min histog. accordingly
+                   //       Bib:
+                   //     2=Marchesini-like (Only Pt sorting: Single sorting)
+                   //        sort only according Pt summation scored, find minimum between both zones and
+                   //        fill Pt and Multiplicity Max and Min summation histog. following only this criterium
+                   //        Bib: Phys. Rev. D 38, 3419 (1988)
+                   //     3=Nameless pt per track single sorting
+                   //        sort according to pt per track scored in each transverse zone 
+                   //        lowest values indicates minimum zone.   
+                   //     4=User Selection sorting (NOTE: USER must implement it within cxx)
+   
     
     // Jet cuts 
-    Double_t   fChJetPtMin;       // Min Pt for charged Particle Jet
-    Double_t   fJet1EtaCut;       // |jet1 eta| < fJet1EtaCut   (fAnaType = 1,2,3)
-    Double_t   fJet2DeltaPhiCut;  // |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut (fAnaType = 2,3)
-    Double_t   fJet2RatioPtCut;   // Jet2.Pt/Jet1Pt > fJet2RatioPtCut  (fAnaType = 2,3)
-    Double_t   fJet3PtCut;        // Jet3.Pt < fJet3PtCut  (fAnaType = 3)
+    Double_t      fChJetPtMin;       // Min Pt for charged Particle Jet
+    Double_t      fJet1EtaCut;       // |jet1 eta| < fJet1EtaCut   (fAnaType = 1,2,3)
+    Double_t      fJet2DeltaPhiCut;  // |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut (fAnaType = 2,3)
+    Double_t      fJet2RatioPtCut;   // Jet2.Pt/Jet1Pt > fJet2RatioPtCut  (fAnaType = 2,3)
+    Double_t      fJet3PtCut;        // Jet3.Pt < fJet3PtCut  (fAnaType = 3)
+
     // track cuts
-    Double_t   fTrackPtCut;       // Pt cut of tracks in the regions
-    Double_t   fTrackEtaCut;      // Eta cut on tracks in the regions (fRegionType=1)
-    Double_t   fAvgTrials;        // average trials used to fill the fh1Triasl histogram in case we do not have trials on a event by event basis
-                
-    // Histograms    ( are owned by fListOfHistos TList )
-    TH1F*  fhNJets;                  //!
-    TH1F*  fhEleadingPt;             //!
-    
-    TH1F*  fhMinRegPtDist;           //!
-    TH1F*  fhRegionMultMin;          //!
-    TH1F*  fhMinRegAvePt;            //!
-    TH1F*  fhMinRegSumPt;            //!
-    TH1F*  fhMinRegMaxPtPart;        //!
-    TH1F*  fhMinRegSumPtvsMult;      //!
-    
-    TH2F*  fhdNdEtaPhiDist;          //!
-    TH2F*  fhFullRegPartPtDistVsEt;  //!
-    TH2F*  fhTransRegPartPtDistVsEt; //!
-    
-    TH1F*  fhRegionSumPtMaxVsEt;     //!
-    TH1I*  fhRegionMultMax;          //!
-    TH1F*  fhRegionMultMaxVsEt;      //!
-    TH1F*  fhRegionSumPtMinVsEt;     //!
-    TH1F*  fhRegionMultMinVsEt;      //!
-    TH1F*  fhRegionAveSumPtVsEt;     //!
-    TH1F*  fhRegionDiffSumPtVsEt;    //!
-    
-    TH1F*  fhRegionAvePartPtMaxVsEt; //!
-    TH1F*  fhRegionAvePartPtMinVsEt; //!
-    TH1F*  fhRegionMaxPartPtMaxVsEt; //!
-    
-    //TH1F*  fhRegForwardSumPtVsEt;    //!
-    //TH1F*  fhRegForwardMultVsEt;     //!
-    //TH1F*  fhRegBackwardSumPtVsEt;   //!
-    //TH1F*  fhRegBackwardMultVsEt;    //!
-    TH2F*  fhRegForwardMult;         //!
-    TH2F*  fhRegForwardSumPtvsMult;  //!
-    TH2F*  fhRegBackwardMult;        //!
-    TH2F*  fhRegBackwardSumPtvsMult; //!
-    TH2F*  fhRegForwardPartPtDistVsEt; //!
-    TH2F*  fhRegBackwardPartPtDistVsEt; //!
-    TH2F*  fhRegTransMult;         //!
-    TH2F*  fhRegTransSumPtVsMult;    //!
-    TH2F*  fhMinRegSumPtJetPtBin;    //!
-    TH2F*  fhMaxRegSumPtJetPtBin;    //!
-    TH1F*  fhVertexMult;
-
-    //        TH2F*  fhValidRegion; //! test to be canceled
-    
-    TProfile* fh1Xsec;               //!
-    TH1F*  fh1Trials;                //!
-    
-    TTree* fSettingsTree;            //! Fast Settings saving
-
-    ClassDef( AliAnalysisTaskUE, 4); // Analysis task for Underlying Event analysis
+    Double_t      fTrackEtaCut;      // Eta cut on tracks in the regions (fRegionType=1)
+    Double_t      fTrackPtCut;       // Pt cut of tracks in the regions
+   
+    // MC cross-section 
+    Double_t      fAvgTrials;        // average trials used to fill the fh1Trials histogram in case we do not have trials on a event by event basis
+    //TProfile*     fh1Xsec;           //!
+    //TH1F*         fh1Trials;         //!
+    
+    ClassDef( AliAnalysisTaskUE, 5); // Analysis task for Underlying Event analysis
   };
 
 #endif