Changes to QA TPC, extra task for cosmics, Added utility classes for UE analysis
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 Jun 2010 10:23:29 +0000 (10:23 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 Jun 2010 10:23:29 +0000 (10:23 +0000)
14 files changed:
PWG4/JetTasks/AliAnalyseUE.cxx [new file with mode: 0644]
PWG4/JetTasks/AliAnalyseUE.h [new file with mode: 0644]
PWG4/JetTasks/AliHistogramsUE.cxx [new file with mode: 0644]
PWG4/JetTasks/AliHistogramsUE.h [new file with mode: 0644]
PWG4/JetTasks/AliPWG4CosmicCandidates.cxx [new file with mode: 0644]
PWG4/JetTasks/AliPWG4CosmicCandidates.h [new file with mode: 0644]
PWG4/JetTasks/AliPWG4HighPtQATPConly.cxx
PWG4/JetTasks/AliPWG4HighPtQATPConly.h
PWG4/JetTasks/AliPWG4HighPtSpectra.cxx
PWG4/JetTasks/AliPWG4HighPtSpectra.h
PWG4/PWG4JetTasksLinkDef.h
PWG4/libPWG4JetTasks.pkg
PWG4/macros/AddTaskPWG4CosmicCandidates.C [new file with mode: 0644]
PWG4/macros/AnalysisTrainPWG4Jets.C

diff --git a/PWG4/JetTasks/AliAnalyseUE.cxx b/PWG4/JetTasks/AliAnalyseUE.cxx
new file mode 100644 (file)
index 0000000..0f22f20
--- /dev/null
@@ -0,0 +1,1205 @@
+/*************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: A.Abrahantes, E.Lopez, S.Vallero                               *
+ * Version 1.0                                                            *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+#include <TROOT.h>
+#include <TBranch.h>
+#include <TCanvas.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH1I.h>
+#include <TH2F.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TMath.h>
+#include <TObjArray.h>
+#include <TProfile.h>
+#include <TRandom.h>
+#include <TSystem.h>
+#include <TTree.h>
+#include <TVector3.h>
+
+#include "AliAnalyseUE.h"
+#include "AliAnalysisTaskUE.h"
+#include "AliAnalysisTask.h"
+#include "AliHistogramsUE.h"
+
+#include "AliAnalysisManager.h"
+#include "AliAODEvent.h"
+#include "AliAODHandler.h"
+#include "AliAODInputHandler.h"
+#include "AliAODJet.h"
+#include "AliAODMCParticle.h"
+#include "AliAODTrack.h"
+#include "AliKFVertex.h"
+#include "AliMCEvent.h"
+#include "AliMCEventHandler.h"
+#include "AliStack.h"
+
+#include "AliAnalysisHelperJetTasks.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliInputEventHandler.h"
+#include "AliLog.h"
+#include "AliStack.h"
+
+////////////////////////////////////////////////
+//--------------------------------------------- 
+// Class for transverse regions analysis
+//---------------------------------------------
+////////////////////////////////////////////////
+
+
+using namespace std;
+
+ClassImp(AliAnalyseUE)
+
+//-------------------------------------------------------------------
+AliAnalyseUE::AliAnalyseUE() :
+  TObject(),
+  //fTaskUE(0),
+  fkAOD(0x0),            
+  fDebug(0),
+  fSimulateChJetPt(kFALSE),
+  fAnaType(1),         
+  fAreaReg(1.5393), // Pi*0.7*0.7
+  fConeRadius(0.7),
+  fFilterBit(0xFF),
+  fRegionType(1),
+  fUseChargeHadrons(kFALSE),
+  fUseChPartJet(kFALSE),
+  fUsePositiveCharge(kTRUE),
+  fUseSingleCharge(kFALSE),
+  fOrdering(1),
+  fJet1EtaCut(0.2),
+  fJet2DeltaPhiCut(2.616),    // 150 degrees
+  fJet2RatioPtCut(0.8),
+  fJet3PtCut(15.),
+  fTrackEtaCut(0.9),
+  fTrackPtCut(0.),
+  fHistos(0x0),
+  fSumPtRegionPosit(0.),
+  fSumPtRegionNegat(0.),
+  fSumPtRegionForward(0.),
+  fSumPtRegionBackward(0.),
+  fMaxPartPtRegion(0.),
+  fNTrackRegionPosit(0),
+  fNTrackRegionNegat(0),
+  fNTrackRegionForward(0),
+  fNTrackRegionBackward(0),
+  fSettingsTree(0x0)
+
+{
+  // constructor
+}
+
+
+//-------------------------------------------------------------------
+AliAnalyseUE::AliAnalyseUE(const AliAnalyseUE & original) :
+  TObject(original),
+  //fTaskUE(original.fTaskUE),
+  fkAOD(original.fkAOD),            
+  fDebug(original.fDebug),
+  fSimulateChJetPt(original.fSimulateChJetPt),
+  fAnaType(original.fAnaType),
+  fAreaReg(original.fAreaReg),
+  fConeRadius(original.fConeRadius),
+  fFilterBit(original.fFilterBit),
+  fRegionType(original.fRegionType),
+  fUseChargeHadrons(original.fUseChargeHadrons),
+  fUseChPartJet(original.fUseChPartJet),
+  fUsePositiveCharge(original.fUsePositiveCharge),
+  fUseSingleCharge(original.fUseSingleCharge),
+  fOrdering(original.fOrdering),
+  fJet1EtaCut(original.fJet1EtaCut),
+  fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),
+  fJet2RatioPtCut(original.fJet2RatioPtCut),
+  fJet3PtCut(original.fJet3PtCut),
+  fTrackEtaCut(original.fTrackEtaCut),
+  fTrackPtCut(original.fTrackPtCut),
+  fHistos(original.fHistos),
+  fSumPtRegionPosit(original.fSumPtRegionPosit),
+  fSumPtRegionNegat(original.fSumPtRegionNegat),
+  fSumPtRegionForward(original.fSumPtRegionForward),
+  fSumPtRegionBackward(original.fSumPtRegionBackward),
+  fMaxPartPtRegion(original.fMaxPartPtRegion),
+  fNTrackRegionPosit(original.fNTrackRegionPosit),
+  fNTrackRegionNegat(original.fNTrackRegionNegat),
+  fNTrackRegionForward(original.fNTrackRegionForward),
+  fNTrackRegionBackward(original.fNTrackRegionBackward),
+  fSettingsTree(original.fSettingsTree)
+{
+  //copy constructor   
+}
+
+//-------------------------------------------------------------------
+AliAnalyseUE & AliAnalyseUE::operator = (const AliAnalyseUE & /*source*/)
+{
+  // assignment operator
+  return *this;
+}
+
+
+//-------------------------------------------------------------------
+AliAnalyseUE::~AliAnalyseUE(){
+
+  //clear memory
+  delete[] fkAOD;
+  fkAOD = NULL;
+
+  
+  
+}
+
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::AnalyseMC(TVector3 *jetVect,AliMCEvent *mcEvent, AliGenPythiaEventHeader *pythiaGenHeader,Int_t conePosition, Bool_t useAliStack, Bool_t constrainDistance, Double_t minDistance){
+
+  // Execute the analysis in case of MC input
+  fSumPtRegionPosit = 0.;
+  fSumPtRegionNegat = 0.;
+  fSumPtRegionForward = 0.;
+  fSumPtRegionBackward = 0.;
+  fMaxPartPtRegion = 0.;
+  fNTrackRegionPosit = 0;
+  fNTrackRegionNegat = 0;
+  fNTrackRegionForward = 0;
+  fNTrackRegionBackward = 0;
+
+  static Double_t const  kPI     = TMath::Pi();
+  static Double_t const  k270rad = 270.*kPI/180.;
+
+  //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 (constrainDistance){
+       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]) > minDistance) 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
+  Double_t maxPtJet1 = random * jetVectnew[index].Pt();  
+    
+  fHistos->FillHistogram("hEleadingPt", maxPtJet1 );    
+    
+  if (useAliStack){//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();
+               Double_t pT = mctrk->Pt();
+               Double_t eta = mctrk->Eta();
+               Int_t pdgCode = mctrk->GetPdgCode();
+
+               if (!TrackMCSelected(charge, pT, eta, pdgCode))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();
+               fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, maxPtJet1 ); 
+        
+               fHistos->FillHistogram("hFullRegPartPtDistVsEt", 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, conePosition );  
+        
+               if (region == 1) {
+                       if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
+                       fSumPtRegionPosit += mctrk->Pt();
+                       fNTrackRegionPosit++;
+                       fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
+                       }
+               if (region == -1) {
+                       if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
+                       fSumPtRegionNegat += mctrk->Pt();
+                       fNTrackRegionNegat++;
+                       fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
+                       }
+               if (region == 2){ //forward
+                       fSumPtRegionForward += mctrk->Pt();
+                       fNTrackRegionForward++;
+                       fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
+                       }
+               if (region == -2){ //backward
+                       fSumPtRegionBackward += mctrk->Pt();
+                       fNTrackRegionBackward++;
+                       fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
+                       }
+               } //end loop on stack particles     
+    }else{//Try mc Particle
+
+      TClonesArray* farray = (TClonesArray*)fkAOD->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;
+        
+       Double_t charge = mctrk->Charge(); 
+       Double_t pT = mctrk->Pt();
+       Double_t eta = mctrk->Eta();
+       Int_t pdgCode = mctrk->GetPdgCode();
+
+       if (!TrackMCSelected(charge, pT, eta, pdgCode))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();
+        fHistos->FillHistogram("hdNdEtaPhiDist", deltaPhi, maxPtJet1 );
+
+       fHistos->FillHistogram("hFullRegPartPtDistVsEt", 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, conePosition );  
+        
+        if (region == 1) { //right
+          if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
+          fSumPtRegionPosit += mctrk->Pt();
+          fNTrackRegionPosit++;
+         fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
+        }
+        if (region == -1) { //left
+          if( fMaxPartPtRegion < mctrk->Pt() ) fMaxPartPtRegion = mctrk->Pt();
+          fSumPtRegionNegat += mctrk->Pt();
+          fNTrackRegionNegat++;
+          fHistos->FillHistogram("hTransRegPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
+        }
+        if (region == 2){ //forward
+          fSumPtRegionForward += mctrk->Pt();
+          fNTrackRegionForward++;
+          fHistos->FillHistogram("hRegForwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
+        }
+        if (region == -2){ //backward
+          fSumPtRegionBackward += mctrk->Pt();
+          fNTrackRegionBackward++;
+         fHistos->FillHistogram("hRegBackwardPartPtDistVsEt", mctrk->Pt(), maxPtJet1 );
+        }
+        
+      }//end loop AliAODMCParticle tracks
+   }  
+}
+
+
+
+//-------------------------------------------------------------------
+Bool_t AliAnalyseUE::AnaTypeSelection(TVector3 *jetVect ){
+
+  // Cut events by jets topology
+  // anaType:
+  //     1 = inclusive,
+  //         - Jet1 |eta| < jet1EtaCut
+  //     2 = back to back inclusive
+  //         - fulfill case 1
+  //         - |Jet1.Phi - Jet2.Phi| > jet2DeltaPhiCut
+  //         - Jet2.Pt/Jet1Pt > jet2RatioPtCut
+  //     3 = back to back exclusive
+  //         - fulfill case 2
+  //         - Jet3.Pt < jet3PtCut
+
+  Double_t eta=jetVect[0].Eta();
+  if( TMath::Abs(eta) > fJet1EtaCut) {
+       if( fDebug > 1 ) AliInfo("\n   Skipping Event...Jet1 |eta| > fJet1EtaCut");
+       return kFALSE;
+       }
+  // back to back inclusive
+  if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() < 0. ) {
+       if( fDebug > 1 ) AliInfo("\n   Skipping Event... no second Jet found");
+       return kFALSE;
+       }
+  if( fAnaType > 1 && fAnaType < 4 && jetVect[1].Pt() > 0. ) {
+       if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
+       jetVect[1].Pt()/jetVect[0].Pt() < fJet2RatioPtCut ) {
+               if( fDebug > 1 ) AliInfo("\n   Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
+               return kFALSE;
+               }
+       }
+  // back to back exclusive
+  if( fAnaType > 2 && fAnaType < 4 && jetVect[2].Pt() > 0. ) {
+       if( jetVect[2].Pt() > fJet3PtCut ) {
+               if( fDebug > 1 ) AliInfo("\n   Skipping Event... Jet3.Pt > fJet3PtCut ");
+               return kFALSE;
+               }
+       }
+  return kTRUE;  
+}
+
+
+//-------------------------------------------------------------------
+TList* AliAnalyseUE::CreateHistograms(Int_t bins, Double_t min, Double_t max){
+  
+  //Initialize histograms from class AliHistogramsUE
+  fHistos = new AliHistogramsUE();
+  TList* list = new TList();
+  list = fHistos->CreateHistos(bins, min, max, fTrackEtaCut);
+
+
+  return list;
+}
+
+
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::FillLeadingJet( Double_t  w){
+
+ fHistos->FillHistogram("hEleadingPt",w);
+
+}
+
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::FillRegions(Bool_t isNorm2Area,  TVector3 *jetVect){
+  
+  // Fill the different topological regions
+  Double_t maxPtJet1 = jetVect[0].Pt();
+  static Double_t const  kPI     = TMath::Pi();
+  static Double_t const  k120rad = 120.*kPI/180.;
+  Double_t const kMyTolerance = 0.0000001;
+
+  //Area for Normalization 
+  // Forward and backward
+  Double_t normArea = 1.;
+  // Transverse
+  if (isNorm2Area) {
+       SetRegionArea(jetVect);
+       normArea =  2.*fTrackEtaCut*k120rad ;
+       } else fAreaReg = 1.;
+  
+  Double_t avePosRegion = (fNTrackRegionPosit) ? fSumPtRegionPosit/fNTrackRegionPosit : 0.;
+  Double_t aveNegRegion = (fNTrackRegionNegat) ? fSumPtRegionNegat/fNTrackRegionNegat : 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( fSumPtRegionPosit > fSumPtRegionNegat ) {
+      FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
+    } else {
+      FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
+    }
+    if (fNTrackRegionPosit > fNTrackRegionNegat ) {
+      FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
+    } else {
+      FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
+    }
+  } else if( fOrdering == 2 ) {
+    if (fSumPtRegionPosit > fSumPtRegionNegat) {
+      FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
+      FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
+    } else {
+      FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
+      FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
+    }
+  } else if( fOrdering == 3 ){
+     if (avePosRegion > aveNegRegion) {
+        FillSumPtRegion( maxPtJet1, fSumPtRegionPosit/fAreaReg, fSumPtRegionNegat/fAreaReg );
+        FillMultRegion( maxPtJet1, fNTrackRegionPosit/fAreaReg, fNTrackRegionNegat/fAreaReg, fSumPtRegionNegat/fAreaReg );
+     }else{
+        FillSumPtRegion( maxPtJet1, fSumPtRegionNegat/fAreaReg, fSumPtRegionPosit/fAreaReg );
+        FillMultRegion( maxPtJet1, fNTrackRegionNegat/fAreaReg, fNTrackRegionPosit/fAreaReg, fSumPtRegionPosit/fAreaReg );
+     }
+  }
+  fHistos->FillHistogram("hRegionMaxPartPtMaxVsEt",maxPtJet1, fMaxPartPtRegion);  
+  
+  // Compute pedestal like magnitudes
+  fHistos->FillHistogram("hRegionDiffSumPtVsEt",maxPtJet1, (TMath::Abs(fSumPtRegionPosit-fSumPtRegionNegat)/(2.0*fAreaReg))+kMyTolerance);
+  fHistos->FillHistogram("hRegionAveSumPtVsEt", maxPtJet1, (fSumPtRegionPosit+fSumPtRegionNegat)/(2.0*fAreaReg));
+
+  // Transverse as a whole
+  fHistos->FillHistogram("hRegTransMult", maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat, (fNTrackRegionPosit + fNTrackRegionNegat)/(2.0*fAreaReg));
+ fHistos->FillHistogram("hRegTransSumPtVsMult",maxPtJet1, fNTrackRegionPosit + fNTrackRegionNegat , (fSumPtRegionNegat + fSumPtRegionPosit)/(2.0 *fAreaReg));
+
+  // Fill Histograms for Forward and away side w.r.t. leading jet direction
+  // Pt dependence
+  //fHistos->FillHistogram("hRegForwardSumPtVsEt",maxPtJet1, fSumPtRegionForward/normArea );
+  //fHistos->FillHistogram("hRegForwardMultVsEt",maxPtJet1, fNTrackRegionForward/normArea );
+  //fHistos->FillHistogram("hRegBackwardSumPtVsEt",maxPtJet1, fSumPtRegionBackward/normArea );
+  //fHistos->FillHistogram("hRegBackwardMultVsEt",maxPtJet1, fNTrackRegionBackward/normArea);
+  
+  // Multiplicity dependence
+  fHistos->FillHistogram("hRegForwardMult", maxPtJet1, fNTrackRegionForward, fNTrackRegionForward/normArea);
+  fHistos->FillHistogram("hRegForwardSumPtvsMult", maxPtJet1, fNTrackRegionForward,fSumPtRegionForward/normArea);
+  fHistos->FillHistogram("hRegBackwardMult", maxPtJet1, fNTrackRegionBackward, fNTrackRegionBackward/normArea );
+  fHistos->FillHistogram("hRegBackwardSumPtvsMult", maxPtJet1, fNTrackRegionBackward,fSumPtRegionBackward/normArea);
+}
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::FillTrials(const char *namex, Double_t  w){
+
+ fHistos->GetTrials()->Fill(namex,w);
+
+}
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::FillVertex(Double_t  w){
+
+ fHistos->FillHistogram("hVertexMult",w);
+
+}
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::FillXsec(const char *namex, Double_t  w){
+
+ fHistos->GetXsec()->Fill(namex,w);
+
+}
+
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::FindMaxMinRegions(TVector3 *jetVect, Int_t conePosition){
+  
+  // Identify the different topological zones
+  fSumPtRegionPosit = 0.;
+  fSumPtRegionNegat = 0.;
+  fSumPtRegionForward = 0.;
+  fSumPtRegionBackward = 0.;
+  fMaxPartPtRegion = 0.;
+  fNTrackRegionPosit = 0;
+  fNTrackRegionNegat = 0;
+  fNTrackRegionForward = 0;
+  fNTrackRegionBackward = 0;
+  static Double_t const  kPI     = TMath::Pi();
+  static Double_t const  kTWOPI  = 2.*kPI;
+  static Double_t const  k270rad = 270.*kPI/180.;
+  Double_t const kMyTolerance = 0.0000001;
+
+    Int_t nTracks = fkAOD->GetNTracks();
+    if (fDebug > 1) AliInfo(Form(" ==== AOD tracks = %d \n ",nTracks));
+    
+    for (Int_t ipart=0; ipart<nTracks; ++ipart) {
+      
+    AliAODTrack* part = fkAOD->GetTrack( ipart );
+    if (fDebug > 1) AliInfo(Form(" ==== AOD track = %d pt = %f charge = %d \n ",ipart,part->Pt(),part->Charge()));
+    // track selection
+    if (! TrackSelected(part)) 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 ) fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, jetVect[0].Pt());
+    else if (TMath::Abs(deltaPhi-k270rad) >= kMyTolerance && TMath::Abs(jetVect[0].Eta()-partVect.Eta()) >= kMyTolerance){
+       fHistos->FillHistogram("hdNdEtaPhiDist",deltaPhi, jetVect[0].Pt());
+       isFlagPart = kFALSE;
+       }
+      
+    fHistos->FillHistogram("hFullRegPartPtDistVsEt", part->Pt(), jetVect[0].Pt());  
+    
+    Int_t region = IsTrackInsideRegion( jetVect, &partVect, conePosition );  
+    if (region == 1) {
+       if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
+        fSumPtRegionPosit += part->Pt();
+        fNTrackRegionPosit++;
+       fHistos->FillHistogram("hTransRegPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
+       }
+    if (region == -1) {
+       if( fMaxPartPtRegion < part->Pt() ) fMaxPartPtRegion = part->Pt();
+        fSumPtRegionNegat += part->Pt();
+        fNTrackRegionNegat++;
+       fHistos->FillHistogram("hTransRegPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
+       }
+    if (region == 2){ //forward
+       fSumPtRegionForward += part->Pt();
+        fNTrackRegionForward++;
+       fHistos->FillHistogram("hRegForwardPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
+       }
+    if (region == -2){ //backward
+       fSumPtRegionBackward += part->Pt();
+        fNTrackRegionBackward++;
+       fHistos->FillHistogram("hRegBackwardPartPtDistVsEt",part->Pt(), jetVect[0].Pt());
+       }
+    }//end loop AOD tracks
+
+}
+
+
+//-------------------------------------------------------------------
+TList* AliAnalyseUE::GetHistograms(){
+
+  TList *list = fHistos->GetListOfHistos();
+
+  return list;
+
+}
+
+
+//-------------------------------------------------------------------
+TVector3 AliAnalyseUE::GetOrderedClusters(TString aodBranch, Bool_t chargedJets, Double_t chJetPtMin){ 
+
+  // jets from AOD, on-the-fly or leading particle
+  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];
+  
+  jetVect[0].SetPtEtaPhi(-1.,-1.,-1.);
+  jetVect[1].SetPtEtaPhi(-1.,-1.,-1.);
+  jetVect[2].SetPtEtaPhi(-1.,-1.,-1.);
+  
+  Int_t nJets = 0;
+  //TClonesArray* fArrayJets;
+  TObjArray* arrayJets;
+  // 1) JETS FROM AOD BRANCH (standard, non-standard or delta)
+  if (!chargedJets && fAnaType != 4 ) { 
+       AliInfo(" ==== Read AODs  !");
+       AliInfo(Form(" ====  Reading Branch: %s  ", aodBranch.Data()));
+       arrayJets = (TObjArray*)fkAOD->GetList()->FindObject(aodBranch.Data());
+       if (!arrayJets){
+              AliFatal(" No jet-array! ");
+              return *jetVect;
+              }
+
+       // Find Leading Jets 1,2,3 
+       // (could be skipped if Jets are sort by Pt...)
+       nJets=arrayJets->GetEntries();
+       for( Int_t i=0; i<nJets; ++i ) {
+               AliAODJet* jet = (AliAODJet*)arrayJets->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*) arrayJets->At(index1);
+               if(jet)jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
+               }
+       if( index2 != -1 ) {
+               AliAODJet* jet= (AliAODJet*) arrayJets->At(index2);
+               if(jet)jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
+               }
+       if( index3 != -1 ) {
+                       AliAODJet* jet = (AliAODJet*) arrayJets->At(index3);
+               if(jet)jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
+               }
+    
+  }
+
+
+  // 2) ON-THE-FLY CDF ALGORITHM
+  if (chargedJets){ 
+    // Printf(" ==== Run CDF algorithm on the fly  !");
+       arrayJets = FindChargedParticleJets(chJetPtMin);
+        if( arrayJets ) {
+               nJets = arrayJets->GetEntriesFast();
+               if( nJets > 0 ) {
+                       index1 = 0;
+                       AliAODJet* jet = (AliAODJet*)arrayJets->At(0);
+                       maxPtJet1 = jet->Pt();
+                       jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
+                       }
+               if( nJets > 1 ) {
+                       index2 = 1;
+                       AliAODJet* jet = (AliAODJet*)arrayJets->At(1);
+                       maxPtJet2 = jet->Pt();
+                       jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
+                       }
+               if( nJets > 2 ) {
+                       index3 = 2;
+                       AliAODJet* jet = (AliAODJet*)arrayJets->At(2);
+                       maxPtJet3 = jet->Pt();
+                       jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz());
+                       }
+      
+               arrayJets->Delete();
+               delete arrayJets;
+               }
+       }
+  
+
+  // 3) LEADING PARTICLE
+  if( fAnaType == 4 ){
+       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; 
+               }
+
+       }
+  //fHistos->FillHistoNJets(nJets);
+  if (fHistos ) fHistos->FillHistogram("hNJets",nJets);
+   
+  return *jetVect;
+
+}
+
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::Initialize(AliAnalysisTaskUE& taskUE){
+   
+  //Get principal settings from current instance of UE analysis-task
+  fAnaType = taskUE.GetAnaTopology();         
+  fkAOD = taskUE.GetAOD();           
+  fConeRadius = taskUE.GetConeRadius();
+  fDebug = taskUE.GetDebugLevel();
+  fFilterBit = taskUE.GetFilterBit();
+  fJet1EtaCut = taskUE.GetJet1EtaCut();
+  fJet2DeltaPhiCut = taskUE.GetJet2DeltaPhiCut();
+  fJet2RatioPtCut = taskUE.GetJet2RatioPtCut();
+  fJet3PtCut = taskUE.GetJet3PtCut();
+  fOrdering = taskUE.GetPtSumOrdering() ;
+  fRegionType = taskUE.GetRegionType();
+  fSimulateChJetPt = taskUE.GetSimulateChJetPt();
+  fTrackEtaCut = taskUE.GetTrackEtaCut(); 
+  fTrackPtCut = taskUE.GetTrackPtCut();
+  fUseChargeHadrons = taskUE.GetUseChargeHadrons();
+  fUseChPartJet = taskUE.GetUseChPartJet();
+  fUsePositiveCharge = taskUE.GetUseNegativeChargeType();
+  fUseSingleCharge = taskUE.GetUseSingleCharge();
+  
+  //Write settings to output list
+  fSettingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
+  fSettingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/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();
+  (fHistos->GetListOfHistos())->Add(fSettingsTree);
+}
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::Initialize(Int_t anaType,AliAODEvent* aod,Double_t coneRadius, Int_t debug, Int_t filterBit, Double_t jet1EtaCut, Double_t jet2DeltaPhiCut, Double_t jet2RatioPtCut, Double_t jet3PtCut, Int_t ordering, Int_t regionType,Bool_t simulateChJetPt, Double_t trackEtaCut, Double_t trackPtCut, Bool_t useChargeHadrons, Bool_t useChPartJet, Bool_t useNegativeChargeType, Bool_t useSingleCharge ){
+   
+  //Get principal settings from generic analysis-task
+  fAnaType = anaType;         
+  fkAOD = aod;           
+  fConeRadius = coneRadius;
+  fDebug = debug;
+  fFilterBit = filterBit;
+  fJet1EtaCut = jet1EtaCut;
+  fJet2DeltaPhiCut = jet2DeltaPhiCut;
+  fJet2RatioPtCut = jet2RatioPtCut;
+  fJet3PtCut = jet3PtCut;
+  fOrdering = ordering ;
+  fRegionType = regionType;
+  fSimulateChJetPt = simulateChJetPt;
+  fTrackEtaCut = trackEtaCut; 
+  fTrackPtCut = trackPtCut;
+  fUseChargeHadrons = useChargeHadrons;
+  fUseChPartJet = useChPartJet;
+  fUsePositiveCharge = useNegativeChargeType;
+  fUseSingleCharge = useSingleCharge;
+  
+  //Write settings to output list
+  fSettingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
+  fSettingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/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();
+  (fHistos->GetListOfHistos())->Add(fSettingsTree);
+}
+
+Bool_t  AliAnalyseUE::VertexSelection(AliAODEvent *aod, Int_t tracks, Double_t zed ){
+
+  //Require 1 vertex (no TPC stand-alone) with a minimum number of tracks and z-coordinate in a limited range
+  Int_t nVertex = aod->GetNumberOfVertices();
+  if( nVertex > 0 ) { // Only one vertex (reject pileup)
+       AliAODVertex* vertex = (AliAODVertex*)aod->GetPrimaryVertex();
+       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 < tracks || TMath::Abs(zVertex) > zed ) {
+               if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
+               return kFALSE;
+               }
+       if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
+       } else {
+               if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
+               return kFALSE;
+               }
+
+  return kTRUE;
+}
+
+// PRIVATE METHODS **************************************************
+
+TObjArray*  AliAnalyseUE::FindChargedParticleJets( Double_t chJetPtMin )
+{
+  // 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 = fkAOD->GetNTracks();
+  if( !nTracks ) return 0;
+  TObjArray tracks(nTracks);
+  
+  for (Int_t ipart=0; ipart<nTracks; ++ipart) {
+       AliAODTrack* part = fkAOD->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() < chJetPtMin) 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;
+}
+
+
+//____________________________________________________________________
+void AliAnalyseUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
+{
+
+  // Fill average particle Pt of control regions
+  
+  // Max cone
+  fHistos->FillHistogram("hRegionAvePartPtMaxVsEt", leadingE, ptMax);
+  // Min cone
+  fHistos->FillHistogram("hRegionAvePartPtMinVsEt", leadingE, ptMin);
+  // MAke distributions for UE comparison with MB data
+  fHistos->FillHistogram("hMinRegAvePt", ptMin);
+
+}
+
+//____________________________________________________________________
+void AliAnalyseUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
+{
+
+  // Fill Nch multiplicity of control regions
+  
+  // Max cone
+  fHistos->FillHistogram("hRegionMultMaxVsEt", leadingE, nTrackPtmax);
+  fHistos->FillHistogram("hRegionMultMax",  nTrackPtmax);
+  
+  // Min cone
+  fHistos->FillHistogram("hRegionMultMinVsEt", leadingE, nTrackPtmin );
+  fHistos->FillHistogram("hRegionMultMin", nTrackPtmin);
+  
+  // MAke distributions for UE comparison with MB data
+  fHistos->FillHistogram("hMinRegSumPtvsMult", nTrackPtmin,ptMin);
+
+}
+
+//____________________________________________________________________
+void AliAnalyseUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
+{
+  // Fill sumPt of control regions
+  
+  // Max cone
+  fHistos->FillHistogram("hRegionSumPtMaxVsEt", leadingE, ptMax);
+  
+  // Min cone
+  fHistos->FillHistogram("hRegionSumPtMinVsEt", leadingE, ptMin);
+  
+  // MAke distributions for UE comparison with MB data
+  fHistos->FillHistogram("hMinRegSumPt", ptMin);
+  fHistos->FillHistogram("hMinRegSumPtJetPtBin", leadingE, ptMin);
+  fHistos->FillHistogram("hMaxRegSumPtJetPtBin", leadingE, ptMax);
+
+}
+
+//-------------------------------------------------------------------
+Int_t AliAnalyseUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect, Int_t conePosition) 
+{  
+  // 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 (conePosition==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 (conePosition==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 (conePosition==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");
+        
+       return region;
+  }
+
+
+
+//-------------------------------------------------------------------
+void  AliAnalyseUE::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);
+    }
+  }
+}
+
+
+
+
+//-------------------------------------------------------------------
+void AliAnalyseUE::SetRegionArea(TVector3 *jetVect)
+{
+  // Set region area
+  Double_t areaCorrFactor=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{
+               areaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) -
+               (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius));
+               fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-areaCorrFactor;
+               }
+       }else AliWarning("Unknown Region 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,areaCorrFactor));
+}
+
+
+//____________________________________________________________________
+TObjArray*  AliAnalyseUE::SortChargedParticles()
+{
+  //  return an array with all charged particles ordered according to their pT .
+  Int_t nTracks = fkAOD->GetNTracks();
+  if( !nTracks ) return 0;
+  TObjArray* tracks = new TObjArray(nTracks);
+
+  for (Int_t ipart=0; ipart<nTracks; ++ipart) {
+       AliAODTrack* part = fkAOD->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;
+  }
+
+
+//____________________________________________________________________
+const Bool_t AliAnalyseUE::TrackMCSelected(Double_t charge, Double_t pT, Double_t eta, Int_t pdgCode){
+  
+  // MC track selection 
+  Double_t const kMyTolerance = 0.0000001;
+  //if (charge == 0. || charge == -99.) return kFALSE;
+  if (charge < kMyTolerance || charge + 99 < kMyTolerance) return kFALSE;
+        
+  if ( fUseSingleCharge ) { // Charge selection
+       if ( fUsePositiveCharge && charge < 0.) return kFALSE; // keep Positives
+       if ( !fUsePositiveCharge && charge > 0.) return kFALSE; // keep Negatives
+       }
+        
+  //Kinematics cuts on particle
+  if ((pT < fTrackPtCut) || (TMath::Abs(eta) > fTrackEtaCut )) return kFALSE;
+        
+  Bool_t isHadron = TMath::Abs(pdgCode)==211 ||
+       TMath::Abs(pdgCode)==2212 ||
+       TMath::Abs(pdgCode)==321;
+        
+  if ( fUseChargeHadrons && !isHadron ) return kFALSE;
+        
+  return kTRUE;
+
+}
+
+//____________________________________________________________________
+const Bool_t AliAnalyseUE::TrackSelected(AliAODTrack* part){
+
+  // Real track selection
+  if ( !part->TestFilterBit(fFilterBit) ) return kFALSE; // track cut selection
+  if ( !part->IsPrimaryCandidate()) return kFALSE; // 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 ) return kFALSE;
+      
+  if ( !part->Charge() ) return kFALSE; //Only charged
+  if ( fUseSingleCharge ) { // Charge selection
+       if ( fUsePositiveCharge && part->Charge() < 0.) return kFALSE; // keep Positives
+       if ( !fUsePositiveCharge && part->Charge() > 0.) return kFALSE; // keep Negatives
+       } 
+    
+  if ( part->Pt() < fTrackPtCut ) return kFALSE;
+  if( TMath::Abs(part->Eta()) > fTrackEtaCut ) return kFALSE;
+
+  return kTRUE;
+}
+
+
+//____________________________________________________________________
+void AliAnalyseUE::WriteSettings(){
+
+  // Print analysis settings
+  if (fDebug>5){
+       AliInfo("All Analysis Settings in Saved Tree");
+       fSettingsTree->Scan();
+       }
+
+}
diff --git a/PWG4/JetTasks/AliAnalyseUE.h b/PWG4/JetTasks/AliAnalyseUE.h
new file mode 100644 (file)
index 0000000..c1f6a70
--- /dev/null
@@ -0,0 +1,150 @@
+//-*- Mode: C++ -*-
+#ifndef ALIANALYSEUE_H
+#define ALIANALYSEUE_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+
+////////////////////////////////////////////////
+//--------------------------------------------- 
+// Class  for transverse regions analysis
+//---------------------------------------------
+////////////////////////////////////////////////
+
+// --- ROOT system ---
+#include <TObject.h> 
+
+#if __GNUC__ >= 3
+using namespace std;
+#endif
+
+class AliAnalysisTaskUE;
+class AliAODEvent;
+class AliAODTrack;
+class AliGenPythiaEventHeader;
+class AliHistogramsUE;
+class AliMCEvent;
+class TVector3;
+
+class AliAnalyseUE : public TObject {
+
+ public: 
+
+  AliAnalyseUE();                                         //constructor
+  AliAnalyseUE(const AliAnalyseUE & g);                   //copy constructor
+  AliAnalyseUE & operator = (const AliAnalyseUE & g);     //assignment operator
+  virtual ~AliAnalyseUE();                                //virtual destructor
+
+  void                 AnalyseMC(TVector3 *jetVect, AliMCEvent *mcEvent, AliGenPythiaEventHeader  *pythiaGenHeader, Int_t conePosition, Bool_t useAliStack, Bool_t constrainDistance, Double_t minDistance);
+  Bool_t       AnaTypeSelection(TVector3 *jetVect);
+
+  TList*       CreateHistograms(Int_t bins, Double_t min, Double_t max);
+
+  void          FillLeadingJet( Double_t  w);
+  void          FillRegions(Bool_t isNorm2Area, TVector3 *jetVect);
+
+  void          FillTrials(const char *namex, Double_t  w);
+  
+  void          FillVertex(Double_t  w);
+  
+  void          FillXsec(const char *namex, Double_t  w);
+
+  void                 FindMaxMinRegions(TVector3 *jetVect, Int_t conePosition);
+  
+  TList*        GetHistograms();
+  
+  TVector3     GetOrderedClusters(TString aodBranch, Bool_t chargedJets, Double_t chJetPtMin);
+
+  void          Initialize(AliAnalysisTaskUE& tmp);
+  void          Initialize(Int_t anaType, AliAODEvent* aod,Double_t coneRadius, Int_t debug, Int_t filterBit, Double_t 
+  jet1EtaCut, Double_t jet2DeltaPhiCut, Double_t jet2RatioPtCut, Double_t jet3PtCut, Int_t ordering, Int_t regionType,Bool_t simulateChJetPt, Double_t trackEtaCut, Double_t trackPtCut, Bool_t useChargeHadrons, Bool_t useChPartJet, Bool_t usePositiveCharge, Bool_t useSingleCharge);
+
+  Bool_t        VertexSelection(AliAODEvent *value, Int_t tracks, Double_t zed);
+  void         WriteSettings();
+
+  // Various setters when you do not want to initialize members from AliAnalysisTaskUE
+  void          SetAnaTopology(Int_t value)            { fAnaType = value; }
+  void          SetAOD(AliAODEvent *value)             { fkAOD = value; }
+  void          SetConeRadius(Double_t value)          { fConeRadius = value; }
+  void          SetDebug(Int_t value)                  { fDebug = value; }
+  void          SetFilterBit(Int_t value)              { fFilterBit = value; }
+  void          SetJet1EtaCut(Double_t value)          { fJet1EtaCut = value; }
+  void          SetJet2DeltaPhiCut(Double_t value)     { fJet2DeltaPhiCut = value; }
+  void          SetJet2RatioPtCut(Double_t value)      { fJet2RatioPtCut = value; }
+  void          SetJet3PtCut(Double_t value)           { fJet3PtCut = value; }
+  void          SetOrdering(Int_t value)               { fOrdering = value; }
+  void          SetRegionType(Int_t value)             { fRegionType = value; }
+  void          SetSimulateChJetPt(Bool_t value)       { fSimulateChJetPt = value; }
+  void          SetTrackEtaCut(Double_t value)         { fTrackEtaCut = value; }
+  void          SetTrackPtCut(Double_t value)          { fTrackPtCut = value; }
+  void          SetUseChargeHadrons(Bool_t value)      { fUseChargeHadrons = value; }
+  void          SetUseChPartJet(Bool_t value)          { fUseChPartJet = value; }
+  void          SetUsePositiveCharge(Bool_t value)     { fUsePositiveCharge = value; }
+  void          SetUseSingleCharge(Bool_t value)       { fUseSingleCharge = value; }
+
+ private:
+
+  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( Double_t chJetPtMin);
+  Int_t        IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect, Int_t conePosition);
+  void         QSortTracks(TObjArray &a, Int_t first, Int_t last);
+  void         SetRegionArea(TVector3 *jetVect);
+  TObjArray*    SortChargedParticles();     
+  virtual const Bool_t TrackSelected(AliAODTrack* part);
+  virtual const Bool_t TrackMCSelected(Double_t charge, Double_t pT, Double_t eta, Int_t pdgCode);
+
+  
+    //AliAnalysisTaskUE    fTaskUE;        //  current instance of the analysis-task
+    const AliAODEvent*   fkAOD;             //! AOD Event 
+    Int_t          fDebug;           //  Debug flag
+
+    
+    // For MC
+    Bool_t         fSimulateChJetPt;      // Naive simulation of charged jet Pt from original Jet in MC Header
+    
+    // Cuts UE analysis
+    Int_t          fAnaType;              // Analysis type on jet topology: 
+    Double_t       fAreaReg;              // Area of the region To be used as normalization factor when filling histograms
+    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)
+    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()
+
+    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:
+    
+    // Jet cuts 
+    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      fTrackEtaCut;      // Eta cut on tracks in the regions (fRegionType=1)
+    Double_t      fTrackPtCut;       // Pt cut of tracks in the regions
+  
+    AliHistogramsUE* fHistos;       // Pointer to histogram class      
+
+    //to fill the different regions
+    Double_t     fSumPtRegionPosit;    // Sum pT in positive region
+    Double_t      fSumPtRegionNegat;   // Sum pT in negative region
+    Double_t      fSumPtRegionForward; // Sum pT in forward region
+    Double_t      fSumPtRegionBackward;        // Sum pT in backward region
+    Double_t      fMaxPartPtRegion;    // Max part pt in region
+    Int_t         fNTrackRegionPosit;  // Tracks in positive region
+    Int_t         fNTrackRegionNegat;  // Tracks in negative region 
+    Int_t         fNTrackRegionForward;        // Track in forward region
+    Int_t         fNTrackRegionBackward;// Tracks in backward region
+
+    //Store analysis settings
+    TTree*       fSettingsTree;        // To store analysis settings
+    
+    ClassDef(AliAnalyseUE,1)
+};
+#endif
diff --git a/PWG4/JetTasks/AliHistogramsUE.cxx b/PWG4/JetTasks/AliHistogramsUE.cxx
new file mode 100644 (file)
index 0000000..307df56
--- /dev/null
@@ -0,0 +1,1052 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: A.Abrahantes, E.Lopez, S.Vallero                               *
+ * Version 1.0                                                            *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id:$ */
+////////////////////////////////////////////////
+//--------------------------------------------- 
+// Class  to handle histograms for UE analysis
+//---------------------------------------------
+////////////////////////////////////////////////
+#include <TROOT.h>
+#include <TBranch.h>
+#include <TCanvas.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH1I.h>
+#include <TH2F.h>
+#include <TLegend.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TMath.h>
+#include <TProfile.h>
+#include <TRandom.h>
+#include <TStyle.h>
+#include <TSystem.h>
+#include <TTree.h>
+#include <TVector3.h>
+
+#include "AliHistogramsUE.h"
+#include "AliAnalyseUE.h"
+#include "AliLog.h"
+
+ClassImp( AliHistogramsUE)
+
+//____________________________________________________________________
+AliHistogramsUE:: AliHistogramsUE():
+TObject(),
+fBinsPtInHist(0),
+fMinJetPtInHist(0.),
+fMaxJetPtInHist(0.),
+fTrackEtaCut(0.),
+fListOfHistos(0x0),  
+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),
+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)
+{
+  // Default constructor
+
+}
+
+//____________________________________________________________________
+AliHistogramsUE:: AliHistogramsUE(TList *list):
+TObject(),
+fBinsPtInHist(0),
+fMinJetPtInHist(0.),
+fMaxJetPtInHist(0.),
+fTrackEtaCut(0.),
+fListOfHistos(0x0),  
+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),
+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)
+{
+  // Constructor, initialize members from list
+       fhNJets = (TH1F*)list->FindObject("hNJets"); 
+       fhEleadingPt = (TH1F*)list->FindObject("hEleadingPt");
+       fhMinRegPtDist = (TH1F*)list->FindObject("hMinRegPtDist");
+       fhRegionMultMin = (TH1F*)list->FindObject("hRegionMultMin");
+       fhMinRegAvePt = (TH1F*)list->FindObject("hMinRegAvePt");
+       fhMinRegSumPt = (TH1F*)list->FindObject("hMinRegSumPt");
+       fhMinRegMaxPtPart = (TH1F*)list->FindObject("hMinRegMaxPtPart");
+       fhMinRegSumPtvsMult = (TH1F*)list->FindObject("hMinRegSumPtvsMult");
+       fhdNdEtaPhiDist = (TH2F*)list->FindObject("hdNdEtaPhiDist");  
+       fhFullRegPartPtDistVsEt = (TH2F*)list->FindObject("hFullRegPartPtDistVsEt");
+       fhTransRegPartPtDistVsEt = (TH2F*)list->FindObject("hTransRegPartPtDistVsEt");
+       fhRegionSumPtMaxVsEt = (TH1F*)list->FindObject("hRegionSumPtMaxVsEt"); 
+       fhRegionMultMax = (TH1I*)list->FindObject("hRegionMultMax");
+       fhRegionMultMaxVsEt = (TH1F*)list->FindObject("hRegionMultMaxVsEt");
+       fhRegionSumPtMinVsEt = (TH1F*)list->FindObject("hRegionSumPtMinVsEt");
+       fhRegionMultMinVsEt = (TH1F*)list->FindObject("hRegionMultMinVsEt");
+       fhRegionAveSumPtVsEt = (TH1F*)list->FindObject("hRegionAveSumPtVsEt");
+       fhRegionDiffSumPtVsEt = (TH1F*)list->FindObject("hRegionDiffSumPtVsEt");
+       fhRegionAvePartPtMaxVsEt = (TH1F*)list->FindObject("hRegionAvePartPtMaxVsEt");
+       fhRegionAvePartPtMinVsEt = (TH1F*)list->FindObject("hRegionAvePartPtMinVsEt");
+       fhRegionMaxPartPtMaxVsEt = (TH1F*)list->FindObject("hRegionMaxPartPtMaxVsEt");
+       fhRegForwardMult = (TH2F*)list->FindObject("hRegForwardMult"); 
+       fhRegForwardSumPtvsMult = (TH2F*)list->FindObject("hRegForwardSumPtvsMult");
+       fhRegBackwardMult = (TH2F*)list->FindObject("hRegBackwardMult");
+       fhRegBackwardSumPtvsMult = (TH2F*)list->FindObject("hRegBackwardSumPtvsMult");
+       fhRegForwardPartPtDistVsEt = (TH2F*)list->FindObject("hRegForwardPartPtDistVsEt");
+       fhRegBackwardPartPtDistVsEt = (TH2F*)list->FindObject("hRegBackwardPartPtDistVsEt");
+       fhRegTransMult = (TH2F*)list->FindObject("hRegTransMult");
+       fhRegTransSumPtVsMult = (TH2F*)list->FindObject("hRegTransSumPtVsMult");
+       fhMinRegSumPtJetPtBin = (TH2F*)list->FindObject("hMinRegSumPtJetPtBin");
+       fhMaxRegSumPtJetPtBin = (TH2F*)list->FindObject("hMaxRegSumPtJetPtBin");
+       fhVertexMult = (TH1F*)list->FindObject("hVertexMult");
+       fh1Xsec = (TProfile*)list->FindObject("h1Xsec");
+       fh1Trials = (TH1F*)list->FindObject("h1Trials"); 
+
+}
+//____________________________________________________________________
+AliHistogramsUE:: AliHistogramsUE(const AliHistogramsUE & original):
+TObject(),
+fBinsPtInHist(original.fBinsPtInHist),
+fMinJetPtInHist(original.fMinJetPtInHist),
+fMaxJetPtInHist(original.fMaxJetPtInHist),
+fTrackEtaCut(original.fTrackEtaCut),
+fListOfHistos(original.fListOfHistos),  
+fhNJets(original.fhNJets),
+fhEleadingPt(original.fhEleadingPt),
+fhMinRegPtDist(original.fhMinRegPtDist),
+fhRegionMultMin(original.fhRegionMultMin),
+fhMinRegAvePt(original.fhMinRegAvePt), 
+fhMinRegSumPt(original.fhMinRegSumPt),            
+fhMinRegMaxPtPart(original.fhMinRegMaxPtPart),
+fhMinRegSumPtvsMult(original.fhMinRegSumPtvsMult),
+fhdNdEtaPhiDist(original.fhdNdEtaPhiDist),        
+fhFullRegPartPtDistVsEt(original.fhFullRegPartPtDistVsEt), 
+fhTransRegPartPtDistVsEt(original.fhTransRegPartPtDistVsEt),
+fhRegionSumPtMaxVsEt(original.fhRegionSumPtMaxVsEt),
+fhRegionMultMax(original.fhRegionMultMax),         
+fhRegionMultMaxVsEt(original.fhRegionMultMaxVsEt),     
+fhRegionSumPtMinVsEt(original.fhRegionSumPtMinVsEt),         
+fhRegionMultMinVsEt(original.fhRegionMultMinVsEt),
+fhRegionAveSumPtVsEt(original.fhRegionAveSumPtVsEt),
+fhRegionDiffSumPtVsEt(original.fhRegionDiffSumPtVsEt),
+fhRegionAvePartPtMaxVsEt(original.fhRegionAvePartPtMaxVsEt),
+fhRegionAvePartPtMinVsEt(original.fhRegionAvePartPtMinVsEt),
+fhRegionMaxPartPtMaxVsEt(original.fhRegionMaxPartPtMaxVsEt),
+fhRegForwardMult(original.fhRegForwardMult),
+fhRegForwardSumPtvsMult(original.fhRegForwardSumPtvsMult),
+fhRegBackwardMult(original.fhRegBackwardMult),
+fhRegBackwardSumPtvsMult(original.fhRegBackwardSumPtvsMult),
+fhRegForwardPartPtDistVsEt(original.fhRegForwardPartPtDistVsEt),
+fhRegBackwardPartPtDistVsEt(original.fhRegBackwardPartPtDistVsEt),
+fhRegTransMult(original.fhRegTransMult),
+fhRegTransSumPtVsMult(original.fhRegTransSumPtVsMult),
+fhMinRegSumPtJetPtBin(original.fhMinRegSumPtJetPtBin),
+fhMaxRegSumPtJetPtBin(original.fhMaxRegSumPtJetPtBin),
+fhVertexMult(original.fhVertexMult),
+fh1Xsec(original.fh1Xsec),
+fh1Trials(original.fh1Trials)
+{
+
+  // Copy constructor
+
+}
+
+
+//______________________________________________________________
+AliHistogramsUE & AliHistogramsUE::operator = (const AliHistogramsUE & /*source*/)
+{
+
+  // assignment operator
+  return *this;
+
+}
+
+//______________________________________________________________
+TObjArray* AliHistogramsUE::CreateCanvas(const Int_t ncanv){
+       
+       // Create canvas for plotting 
+        printf("Creating %d canvas ... \n",ncanv);
+        TObjArray *arr=new TObjArray;
+       TString name;
+        for(Int_t i=0;i<ncanv;i++ ){
+        name=Form("Canvas %d",i);
+        TCanvas *c=new TCanvas(name,name);
+        c->SetFillColor(0);
+       gStyle->SetOptStat(0);
+       gStyle->SetOptTitle(0);
+       arr->Add(c);
+        }
+return arr;
+}
+
+//____________________________________________________________________
+TList*  AliHistogramsUE::CreateHistos(Int_t bins, Double_t min, Double_t max, Double_t etacut)
+{
+
+  // Create all histograms necessary for UE analysis 
+  fBinsPtInHist = bins;
+  fMinJetPtInHist = min;
+  fMaxJetPtInHist = max;
+  fTrackEtaCut= etacut;
+
+  fListOfHistos = new TList();
+
+  //Number of reconstructed clusters  
+  fhNJets = new TH1F("hNJets", "Number of clusters",  20, 0, 20);
+  fhNJets->SetXTitle("Number of reconstructed clusters");
+  fhNJets->SetYTitle("#");
+  fhNJets->Sumw2();
+  fListOfHistos->Add( fhNJets );                 // At(0) 
+  
+  //pT distribution of leading clusters
+  fhEleadingPt  = new TH1F("hEleadingPt",   "Leading cluster p_{T}",  bins, min,   max);
+  fhEleadingPt->SetXTitle("p_{T} of  cluster (GeV/c)");
+  fhEleadingPt->SetYTitle("1/N_{ev} dN/dp_{T} (|#eta|<0.5)");
+  fhEleadingPt->Sumw2();
+  fListOfHistos->Add( fhEleadingPt );            // At(1)
+  
+  //Track pT distribution in MIN zone
+  fhMinRegPtDist = new TH1F("hMinRegPtDist",   "p_{T} distribution in MIN zone",  50,0.,20.);
+  fhMinRegPtDist->SetXTitle("Track p_{T} (GeV/c)");
+  fhMinRegPtDist->SetYTitle("dN/dp_{T}");
+  fhMinRegPtDist->Sumw2();
+  fListOfHistos->Add( fhMinRegPtDist );          // At(2)
+  
+  //Multiplicity in MIN zone
+  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)            
+  
+  //Average pT in MIN region
+  fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT",  50, 0.,   20.);
+  fhMinRegAvePt->SetXTitle("p_{T} (GeV/c)");
+  fhMinRegAvePt->Sumw2();
+  fListOfHistos->Add( fhMinRegAvePt );           // At(4)
+  
+  //Sum pT in MIN region
+  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)
+  
+  //Track with maximum pT in MIN region
+  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)
+  
+  //Sum pT vs. multiplicity in MIN region
+  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);
+  
+  //Phi correlation track-cluster vs. leading cluster pT 
+  fhdNdEtaPhiDist  = new TH2F("hdNdEtaPhiDist", Form("Charge particle density |#eta|<%3.1f vs #Delta#phi", fTrackEtaCut),62, 0.,   2.*TMath::Pi(), bins, min, max);
+  fhdNdEtaPhiDist->SetXTitle("#Delta#phi");
+  fhdNdEtaPhiDist->SetYTitle("Leading cluster p_{T}");
+  fhdNdEtaPhiDist->Sumw2();
+  fListOfHistos->Add( fhdNdEtaPhiDist );        // At(8)
+  
+  //Can be used to get track pT distribution for different cluster pT bins (full region)
+  fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", Form( "dN/dp_{T} |#eta|<%3.1f vs Leading cluster p_{T}", fTrackEtaCut),100,0.,50., bins, min, max);
+  fhFullRegPartPtDistVsEt->SetYTitle("Leading cluster p_{T}");
+  fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
+  fhFullRegPartPtDistVsEt->Sumw2();
+  fListOfHistos->Add( fhFullRegPartPtDistVsEt );  // At(9) 
+  
+  //Can be used to get part pT distribution for different cluster pT bins (transverse region)
+  fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", Form( "dN/dp_{T} in tranvese regions |#eta|<%3.1f vs Leading cluster p_{T}", fTrackEtaCut),100,0.,50., bins, min,   max);
+  fhTransRegPartPtDistVsEt->SetYTitle("Leading cluster p_{T}");
+  fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
+  fhTransRegPartPtDistVsEt->Sumw2();
+  fListOfHistos->Add( fhTransRegPartPtDistVsEt );  // At(10)
+  
+  //Sum pT in MAX region vs. leading-cluster pT
+  fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt",  "P_{T}^{90, max} vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionSumPtMaxVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegionSumPtMaxVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionSumPtMaxVsEt );      // At(11)
+  
+
+  //Sum pT in MIN region vs. leading-cluster pT
+  fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt",   "P_{T}^{90, min} vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionSumPtMinVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegionSumPtMinVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionSumPtMinVsEt );      // At(12)
+  //Multiplicity in MAX region
+  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)
+  
+  //Multiplicity in MAX region vs. leading-cluster pT
+  fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt",  "N_{ch}^{90, max} vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionMultMaxVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegionMultMaxVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionMultMaxVsEt );       // At(14)
+  
+  //Multiplicity in MIN region vs. leading-cluster pT
+  fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt",  "N_{ch}^{90, min} vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionMultMinVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegionMultMinVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionMultMinVsEt );      // At(15)
+  //Average sum pT in TRANSVERSE(MIN+MAX) region vs. leading-cluster pT 
+  fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionAveSumPtVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegionAveSumPtVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionAveSumPtVsEt );     // At(16)
+  
+  //Difference sum pT (MAX-MIN) vs.  leading-cluster pT
+  fhRegionDiffSumPtVsEt= new TH1F("hRegionDiffSumPtVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
+  fhRegionDiffSumPtVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionDiffSumPtVsEt );    // At(17)
+  
+  //Average track pT in MAX region vs. leading-cluster pT
+  fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionAvePartPtMaxVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegionAvePartPtMaxVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionAvePartPtMaxVsEt );  // At(18)
+  
+  //Average track pT in MIN region vs. leading-cluster pT
+  fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionAvePartPtMinVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegionAvePartPtMinVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionAvePartPtMinVsEt );   // At(19)
+  
+  //Maximum track pT in MAX region vs. leading-cluster pT
+  fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading cluster p_{T}",  bins, min,   max);
+  fhRegionMaxPartPtMaxVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegionMaxPartPtMaxVsEt->Sumw2();
+  fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt );    // At(20)
+  
+  //Multiplicity in FORWARD region
+  fhRegForwardMult = new TH2F("hRegForwardMult", "N_{ch}^{forward}", bins, min, max, 21, -0.5,   20.5);
+  fhRegForwardMult->SetXTitle("N_{ch tracks}");
+  fhRegForwardMult->Sumw2();
+  fListOfHistos->Add( fhRegForwardMult );           // At(25)
+  
+  //Sum pT in FORWARD region vs. multiplicity
+  fhRegForwardSumPtvsMult = new TH2F("hRegForwardSumPtvsMult", "Forward #Sigma p_{T} vs. Multiplicity ", bins, min, max, 21, -0.5,   20.5);
+  fhRegForwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
+  fhRegForwardSumPtvsMult->SetXTitle("N_{charge}");
+  fhRegForwardSumPtvsMult->Sumw2();
+  fListOfHistos->Add( fhRegForwardSumPtvsMult );     // At(26);
+  
+  //Multiplicity in BACKWARD region
+  fhRegBackwardMult = new TH2F("hRegBackwardMult", "N_{ch}^{backward}", bins, min, max, 21, -0.5,   20.5);
+  fhRegBackwardMult->SetXTitle("N_{ch tracks}");
+  fhRegBackwardMult->Sumw2();
+  fListOfHistos->Add( fhRegBackwardMult );           // At(27)
+  //Sum pT in BACKWARD region vs. multiplicity
+  fhRegBackwardSumPtvsMult = new TH2F("hRegBackwardSumPtvsMult", "Backward #Sigma p_{T} vs. Multiplicity ", bins, min, max, 21, -0.5,   20.5);
+  fhRegBackwardSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
+  fhRegBackwardSumPtvsMult->SetXTitle("N_{charge}");
+  fhRegBackwardSumPtvsMult->Sumw2();
+  fListOfHistos->Add( fhRegBackwardSumPtvsMult );     // At(28);
+  
+  //Track pT distribution in FORWARD region vs. leading-cluster pT 
+  fhRegForwardPartPtDistVsEt = new TH2F("hRegForwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading cluster p_{T}", fTrackEtaCut), 100,0.,50., bins, min, max);
+  fhRegForwardPartPtDistVsEt->SetYTitle("Leading cluster p_{T}");
+  fhRegForwardPartPtDistVsEt->SetXTitle("p_{T} (GeV/c)");
+  fhRegForwardPartPtDistVsEt->Sumw2();
+  fListOfHistos->Add( fhRegForwardPartPtDistVsEt );  // At(29) 
+  
+  //Track pT distribution in BACKWARD region vs. leading-cluster pT 
+  fhRegBackwardPartPtDistVsEt = new TH2F("hRegBackwardPartPtDistVsEt", Form( "dN/dP_{T} |#eta|<%3.1f vs Leading cluster p_{T}", fTrackEtaCut), 100,0.,50., bins, min, max);
+  fhRegBackwardPartPtDistVsEt->SetYTitle("Leading cluster p_{T}");
+  fhRegBackwardPartPtDistVsEt->SetXTitle("p_{T}");
+  fhRegBackwardPartPtDistVsEt->Sumw2();
+  fListOfHistos->Add( fhRegBackwardPartPtDistVsEt );  // At(30) 
+  
+  //Multiplicity in TRANSVERSE (MIN+MAX) region
+  fhRegTransMult  = new TH2F("hRegTransMult", "N_{ch}^{transv}", bins, min, max, 21, -0.5,   20.5);
+  fhRegTransMult->SetXTitle("N_{ch tracks}");
+  fhRegTransMult->Sumw2();
+  fListOfHistos->Add( fhRegTransMult );           // At(31)
+  
+  //Sum pT in TRANSVERSE (MIN+MAX) region vs. multiplicity
+  fhRegTransSumPtVsMult = new TH2F("hRegTransSumPtVsMult", "Transverse #Sigma p_{T} vs. Multiplicity ",bins, min, max, 21, -0.5,   20.5);
+  fhRegTransSumPtVsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
+  fhRegTransSumPtVsMult->SetXTitle("N_{charge}");
+  fhRegTransSumPtVsMult->Sumw2();
+  fListOfHistos->Add( fhRegTransSumPtVsMult );     // At(32);
+  
+  //Sum pT in MIN region per cluster pT bin
+  fhMinRegSumPtJetPtBin = new TH2F("hMinRegSumPtJetPtBin", "Transverse Min Reg #Sigma p_{T} per cluster pT bin",bins, min, max, 50, 0.,   20.);
+  fhMinRegSumPtJetPtBin->SetXTitle("Leading cluster p_{T}");
+  fhMinRegSumPtJetPtBin->Sumw2();
+  fListOfHistos->Add( fhMinRegSumPtJetPtBin );           // At(33)
+  
+  //Sum pT in MAX region per cluster pT bin
+  fhMaxRegSumPtJetPtBin = new TH2F("hMaxRegSumPtJetPtBin",      "Transverse Max Reg #Sigma p_{T} per cluster pT bin", bins, min, max, 50, 0.,   20.);
+  fhMaxRegSumPtJetPtBin->SetXTitle("Leading cluster p_{T}");
+  fhMaxRegSumPtJetPtBin->Sumw2();
+  fListOfHistos->Add( fhMaxRegSumPtJetPtBin );           // At(34)
+  
+  //Multiplicity in main vertex
+  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("h1Xsec","xsec from pyxsec.root",1,0,1); 
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+  fh1Xsec->Sumw2();
+  fListOfHistos->Add( fh1Xsec );            //At(36)
+  
+  fh1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+  fh1Trials->Sumw2();
+  fListOfHistos->Add( fh1Trials ); //At(37)
+
+  return fListOfHistos;
+}
+
+
+//____________________________________________________________________
+void AliHistogramsUE::DrawUE(Int_t debug){
+
+    // To draw histograms at the end of task running 
+    // 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!!!!!
+    
+    Int_t binsPtInHist = fhEleadingPt->GetNbinsX();
+    Double_t minJetPtInHist = fhEleadingPt->GetXaxis()->GetBinLowEdge(1);
+    Double_t maxJetPtInHist = fhEleadingPt->GetXaxis()->GetBinUpEdge(binsPtInHist);
+     
+    //Sum pT
+    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" , "", binsPtInHist,  minJetPtInHist, maxJetPtInHist);
+    //TH1F *h1r = new TH1F();
+    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" , "", binsPtInHist,  minJetPtInHist, maxJetPtInHist);
+    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" , "", binsPtInHist,  minJetPtInHist, maxJetPtInHist);
+    //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");
+
+    c1->cd(4);
+    TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading",   "",  binsPtInHist, minJetPtInHist,   maxJetPtInHist);
+    TH1F *h6r = new TH1F("hRegionMultMinVsEtleading",   "",  binsPtInHist, minJetPtInHist,   maxJetPtInHist);
+    h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
+    h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
+    //h5r->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
+    //Double_t xsec = fh1Xsec->GetBinContent(1);
+    Double_t xsec = fh1Xsec->GetBinContent(1);
+    Double_t ntrials = fh1Trials->GetBinContent(1);
+    Double_t normFactor = xsec/ntrials;
+    if(debug > 1)Printf("xSec %f nTrials %f Norm %f \n",xsec,ntrials,normFactor);
+    
+    
+    //Jet pT distribution
+    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(minJetPtInHist);
+    Int_t xbin2 = fhdNdEtaPhiDist->GetYaxis()->FindFixBin(maxJetPtInHist);
+    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);
+    c2->Update(); 
+    
+    //pT distributions
+    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();
+    
+    
+
+    //Save canvas
+    c1->SaveAs("c1.pdf");
+    AliInfo("Canvas 1 saved");
+    c2->SaveAs("c2.pdf");
+    AliInfo("Canvas 2 saved");
+    c3->SaveAs("c3.pdf");
+    AliInfo("Canvas 3 saved");
+  
+}
+
+//____________________________________________________________________
+void AliHistogramsUE::FillHistogram(const char* name, Double_t fillX){
+  
+  // Fill 1D histogram with double
+  ((TH1F*)fListOfHistos->FindObject(name))->Fill(fillX);
+
+}
+
+//____________________________________________________________________
+void AliHistogramsUE::FillHistogram(const char* name, Int_t fillX){
+  
+  // Fill 1D histogram with integer
+  ((TH1F*)fListOfHistos->FindObject(name))->Fill(fillX);
+
+}
+
+//____________________________________________________________________
+void AliHistogramsUE::FillHistogram(const char* name, Double_t fillX, Double_t fillY){
+  
+ // Case of TH1F with weight or TH2F w/o weight
+ TObject *obj = fListOfHistos->FindObject(name);
+ if (obj->InheritsFrom("TH1F")){
+       ((TH1F*)fListOfHistos->FindObject(name))->Fill(fillX, fillY);
+ } else {
+       ((TH2F*)fListOfHistos->FindObject(name))->Fill(fillX, fillY);
+       }
+
+
+}
+
+//____________________________________________________________________
+void AliHistogramsUE::FillHistogram(const char* name, Double_t fillX, Double_t fillY, Double_t weight){
+
+  // Fill 2D histogram with double and weight   
+  ((TH2F*)fListOfHistos->FindObject(name))->Fill(fillX, fillY, weight);
+
+}
+
+//____________________________________________________________________
+void AliHistogramsUE::FillHistogram(const char* name, Double_t fillX, Int_t fillY, Double_t weight){
+  
+  // Fill 2D histogram with integer and weight   
+  ((TH2F*)fListOfHistos->FindObject(name))->Fill(fillX, fillY, weight);
+
+}
+
+//____________________________________________________________________
+TObjArray*  AliHistogramsUE::GetHistosForPlotting(TString data, TString branches){
+
+        // Instance filled histos for plotting purpose
+       printf("Creating histograms ... \n");
+               
+       printf("Reading file: %s\n",data.Data());
+       
+        // Read input files ----------------------------------------- 
+        TFile *fdata = new TFile(data.Data());
+       TDirectoryFile *ddata[20];
+       TList *ldata[20];
+
+       TObjArray  *arrb=branches.Tokenize(";");
+       TIter next(arrb);
+       TObject *o=0;
+       Int_t br=0;
+       while ( (o=next()) ){
+               ddata[br] = (TDirectoryFile*)fdata->Get(Form("PWG4_UE%s",o->GetName()));
+               if(!ddata[br]) printf("ERROR: No histo dir found! \n");
+               ldata[br] = (TList*)ddata[br]->Get(Form("histosUE%s",o->GetName()));
+               printf("Reading branch: %s\n",o->GetName());
+               if(!ldata[br]) printf("ERROR: No histo list found! \n");
+               br++;
+       }
+       
+       TObjArray *arr=new TObjArray;
+
+       TH1F *hjets[20];                // accepted leading jets
+        TH1F *hnjets[20];          // number of accepted jets
+       TH2F *hetaphi[20];         // delta-phi particle-jet correlation
+       TH2F *hptfull[20];         // particle pT all regions vs. jet pT 
+       TH2F *hpttransv[20];       // particle pT transv. regions vs. jet pT 
+       TH1F *hmax[20];         // sum pT in MAX region
+        TH1F *hmin[20];        // sum pT in MIN region
+        TH1F *hmultmax[20];    // multiplicity in MAX region
+        TH1F *hmultmin[20];       // multiplicity in MIN region
+
+       for (Int_t i =0; i<br; i++){
+
+               //Number of jets --------------------------------------------
+               hnjets[i] = (TH1F*) ldata[i]->FindObject("fhNJets"); 
+                       hnjets[i]->GetXaxis()->SetTitle("Number of jets");
+               hnjets[i]->GetYaxis()->SetTitle("1/n_{ev} dN/dN_{jets}");
+               hnjets[i]->SetMarkerStyle(20); 
+               hnjets[i]->SetMarkerColor(1+i);
+
+
+               //Leading jets ----------------------------------------------
+               hjets[i] = (TH1F*) ldata[i]->FindObject("hEleadingPt"); 
+               hjets[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
+               hjets[i]->GetYaxis()->SetTitle("1/n_{ev} dN/dp_{T} (|#eta<0.5|)");
+               hjets[i]->SetMarkerStyle(20); 
+               hjets[i]->SetMarkerColor(1+i);
+               hjets[i]->SetMinimum(0.1);
+               hjets[i]->SetMaximum(1000.);
+
+
+               //Transverse Region MAX -------------------------------------
+               hmax[i] = (TH1F*) ldata[i]->FindObject("hRegionSumPtMaxVsEt");
+               if (!hmax[i])AliInfo("Histo not found!!!");
+               hmax[i]->GetXaxis()->SetTitle("Leading jet P_{T} (GeV/c)");
+               hmax[i]->GetYaxis()->SetTitle("P_{T}^{90,max} (GeV/c)");
+               hmax[i]->SetMarkerStyle(20);
+               hmax[i]->SetMarkerColor(1+i);
+               hmax[i]->Divide(hjets[i]); // normalize for jet spectrum
+               hmax[i]->SetMaximum(5.);
+
+
+               //Transverse Region MIN -------------------------------------
+               hmin[i] = (TH1F*) ldata[i]->FindObject("hRegionSumPtMinVsEt");
+               hmin[i]->GetXaxis()->SetTitle("Leading jet P_{T} (GeV/c)");
+               hmin[i]->GetYaxis()->SetTitle("P_{T}^{90,min} (GeV/c)");
+               hmin[i]->SetMarkerStyle(20);
+               hmin[i]->SetMarkerColor(1+i);
+               hmin[i]->SetMaximum(3.);
+               hmin[i]->Divide(hjets[i]); // normalize for jet spectrum
+
+
+               //Multiplicity MAX ------------------------------------------
+               hmultmax[i] = (TH1F*) ldata[i]->FindObject("hRegionMultMaxVsEt");
+               hmultmax[i]->GetXaxis()->SetTitle("Leading Jet P_{T} (GeV/c)");
+               hmultmax[i]->GetYaxis()->SetTitle("N_{ch}^{90,max}");
+               hmultmax[i]->SetMarkerStyle(20);
+               hmultmax[i]->SetMarkerColor(1+i);
+               hmultmax[i]->SetMaximum(10.);
+               hmultmax[i]->Divide(hjets[i]); // normalize for jet spectrum
+
+
+               //Multiplicity MIN ------------------------------------------
+               hmultmin[i] = (TH1F*) ldata[i]->FindObject("hRegionMultMinVsEt");
+               hmultmin[i]->GetXaxis()->SetTitle("Leading Jet P_{T} (GeV/c)");
+               hmultmin[i]->GetYaxis()->SetTitle("N_{ch}^{90,min}");
+               hmultmin[i]->SetMarkerStyle(20);
+               hmultmin[i]->SetMarkerColor(1+i);
+               hmultmin[i]->SetMaximum(3.);
+               hmultmin[i]->Divide(hjets[i]); // normalize for jet spectrum
+            
+
+               // Phi-correlation with leading jet --------------------------
+               hetaphi[i] = (TH2F*) ldata[i]->FindObject("hdNdEtaPhiDist");
+               hetaphi[i]->GetXaxis()->SetTitle("#Delta #phi (w.r.t. leading jet)");
+               hetaphi[i]->SetMarkerStyle(20);
+            
+
+               // pT distribution in full region vs. jet pT --------------------------
+               hptfull[i] = (TH2F*) ldata[i]->FindObject("hFullRegPartPtDistVsEt");
+               hptfull[i]->GetYaxis()->SetTitle("Leading Jet P_{T} (GeV/c)");
+               hptfull[i]->GetXaxis()->SetTitle("Track P_{T} (GeV/c)");
+            
+
+               // pT distribution in transv region vs. jet pT --------------------------
+               hpttransv[i] = (TH2F*) ldata[i]->FindObject("hTransRegPartPtDistVsEt");
+               hpttransv[i]->GetYaxis()->SetTitle("Leading Jet P_{T} (GeV/c)");
+               hpttransv[i]->GetXaxis()->SetTitle("Track P_{T} (GeV/c)");
+            
+
+               // Return Histos
+               arr->Add(hnjets[i]);      //at 0 number of jets
+               arr->Add(hjets[i]);     //at 1 leading jets
+               arr->Add(hmax[i]);      //at 2 sum pT MAX
+               arr->Add(hmin[i]);      //at 3 sumpT MIN
+               arr->Add(hmultmax[i]);   //at 4 multiplicity MAX 
+               arr->Add(hmultmin[i]);   //at 5 multiplicity MIN
+               arr->Add(hetaphi[i]);     //at 6 phi correlation
+               arr->Add(hptfull[i]);     //at 7 pT distr in full region
+               arr->Add(hpttransv[i]);   //at 8 pT distr in transv region
+
+       }
+       return arr;
+}
+
+//_______________________________________________________________________________________
+void AliHistogramsUE::SetStyle(){
+
+  // Set plotting style
+  gPad->SetFrameFillColor(0);
+  gPad->SetFillColor(0);
+  gPad->SetBorderSize(2);
+  gPad->SetGridy();
+  gPad->SetFrameBorderMode(0);
+  //gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+
+}
+
+
+//____________________________________________________________________
+TList* AliHistogramsUE::GetListOfHistos(){
+
+  // Return list of relevant histograms 
+  return fListOfHistos;
+
+}
+
+//____________________________________________________________________
+void AliHistogramsUE::PlotBranchesUE(TString file, TString branches, Double_t minJetProjection){
+
+  // Function to be called by external macro to plot analysis from different jet branches 
+  
+  //Count the branches 
+  TObjArray  *arr=branches.Tokenize(";");
+  TIter next(arr);
+  TObject *o=0;
+  Int_t br=0;
+  while ( (o=next()) ){
+       br++;
+        }
+  
+
+  //Canvas
+  TObjArray *arrC=CreateCanvas(9);
+  //Histograms 
+  TObjArray *arrH=GetHistosForPlotting(file.Data(),branches.Data());
+  TH1F *hnjets[20];
+  TH1F *hjets[20];
+  TH1F *hmax[20];
+  TH1F *hmin[20];
+  TH1F *hmultmax[20];
+  TH1F *hmultmin[20];
+  TH2F *hetaphi[20];
+  TH2F *hptfull[20];
+  TH2F *hpttransv[20];
+
+  for (Int_t i= 0; i<br; i++){
+       hnjets[i]=((TH1F*)arrH->At(9*i));
+       hjets[i]=((TH1F*)arrH->At(1+(9*i)));
+       hmax[i]=((TH1F*)arrH->At(2+(9*i)));
+       hmin[i]=((TH1F*)arrH->At(3+(9*i)));
+       hmultmax[i]=((TH1F*)arrH->At(4+(9*i)));
+       hmultmin[i]=((TH1F*)arrH->At(5+(9*i)));
+       hetaphi[i]=((TH2F*)arrH->At(6+(9*i))); 
+       hptfull[i]=((TH2F*)arrH->At(7+(9*i))); 
+       hpttransv[i]=((TH2F*)arrH->At(8+(9*i))); 
+       }
+
+  //Define jet-pT range in projections
+  Int_t binstartprojection = hetaphi[0]->GetYaxis()->FindFixBin(minJetProjection); //be careful...  
+  Int_t binstopprojection = hetaphi[0]->GetYaxis()->GetNbins();
+
+  Double_t normphi[20]; 
+
+  for (Int_t i= 0; i<br; i++){
+       normphi[i]=hjets[i]->Integral(binstartprojection,binstopprojection);    
+       hnjets[i]->Scale(1./(hnjets[i]->GetBinWidth(1)*hnjets[i]->GetEntries()));
+        hjets[i]->Scale(1./(hjets[i]->GetBinWidth(1)*hnjets[i]->GetEntries()));
+        }
+
+  //LEGENDS
+  //Legend jets
+  TLegend *leg=new TLegend(0.5,0.6,0.89,0.89);  // for jet spectrum
+  leg->SetFillColor(0);
+  leg->SetHeader("Jet Finders:");
+  //Legend density
+  TLegend *legd = new TLegend(0.1077586,0.6016949,0.4971264,0.8919492,NULL,"brNDC");
+  legd->SetFillColor(0);
+  legd->SetHeader("Jet Finders:");
+  //Legend pT distributions
+  TLegend *legp = new TLegend(0.1364943,0.1292373,0.5258621,0.4194915,NULL,"brNDC");
+  legp->SetFillColor(0);
+  legp->SetHeader("Jet Finders:");
+
+  arr=branches.Tokenize(";");
+  TIter next1(arr);
+  o=0;
+  Int_t brleg=0;
+  while ( (o=next1()) ){
+       leg->AddEntry(hjets[brleg],Form("UE%s",o->GetName()),"p");
+       legd->AddEntry(hjets[brleg],Form("UE%s",o->GetName()),"p");
+       legp->AddEntry(hjets[brleg],Form("UE%s",o->GetName()),"p");
+        brleg++;
+        }
+
+  //1) NUMBER OF CLUSTERS 
+  TCanvas *c0=((TCanvas*)arrC->At(0)); 
+  c0->cd();
+  SetStyle();
+  gPad->SetLogy();
+  hnjets[0]->SetMaximum(1.4);
+  hnjets[0]->Draw("E1");
+  for (Int_t i= 1; i<br; i++){
+       hnjets[i]->Draw("same");
+        }
+  leg->Draw("same");
+
+  //2) LEADING CLUSTERS pT 
+  TCanvas *c1=((TCanvas*)arrC->At(1)); 
+  c1->cd();
+  SetStyle();
+  gPad->SetLogy();
+  hjets[0]->Draw("E1");
+  hjets[0]->SetMaximum(1.4);
+  for (Int_t i= 1; i<br; i++){
+       hjets[i]->Draw("same");
+        }
+  leg->Draw("same");
+
+  //3) SUM-pT IN MAX REGION
+  TCanvas *c2=((TCanvas*)arrC->At(2)); 
+  c2->cd();
+  SetStyle();
+  hmax[0]->Draw("E1");
+  for (Int_t i= 1; i<br; i++){
+       hmax[i]->Draw("same");
+        }
+  leg->Draw("same");
+
+  //4) SUM-pT IN MIN REGION
+  TCanvas *c3=((TCanvas*)arrC->At(3)); 
+  c3->cd();
+  SetStyle();
+  hmin[0]->GetYaxis()->SetRangeUser(0.,2.5);
+  hmin[0]->Draw("E1");
+  for (Int_t i= 1; i<br; i++){
+       hmin[i]->Draw("same");
+        }
+  leg->Draw("same");
+
+  //5) MULTIPLICITY IN MAX REGION
+  TCanvas *c4=((TCanvas*)arrC->At(4));
+  c4->cd();
+  SetStyle();
+  hmultmax[0]->GetYaxis()->SetRangeUser(0.,5.);
+  hmultmax[0]->Draw("E1");
+  for (Int_t i= 1; i<br; i++){
+       hmultmax[i]->Draw("same");
+        }
+  leg->Draw("same");
+
+  //6) MULTIPLICITY IN MIN REGION
+  TCanvas *c5=((TCanvas*)arrC->At(5));
+  c5->cd();
+  SetStyle();
+  hmultmin[0]->GetYaxis()->SetRangeUser(0.,2.5);
+  hmultmin[0]->Draw("E1");
+  for (Int_t i= 1; i<br; i++){
+       hmultmin[i]->Draw("same");
+        }
+  leg->Draw("same");
+
+  //7) JET-TRACK CORRELATION
+  TCanvas *c6=((TCanvas*)arrC->At(6));
+  c6->cd();
+  SetStyle();
+  gPad->SetLogy();
+  TH1D *tmpetaphi[20];
+  for (Int_t i= 0; i<br; i++){
+       tmpetaphi[i]=new TH1D();
+       tmpetaphi[i]=hetaphi[i]->ProjectionX(Form("data%d",i),binstartprojection);
+       tmpetaphi[i]->GetYaxis()->SetTitle("1/n_{lj} dN/d#Delta #phi");
+       tmpetaphi[i]->Scale(1./(hetaphi[i]->GetBinWidth(1)*normphi[i]));
+       tmpetaphi[i]->SetMarkerColor(i+1);
+       tmpetaphi[i]->GetXaxis()->SetLimits(-3.*TMath::Pi()/2.,TMath::Pi()/2.);
+       tmpetaphi[i]->GetYaxis()->SetRangeUser(0.5,20.);
+       if (i==0) tmpetaphi[i]->Draw("E1");
+       else tmpetaphi[i]->Draw("same");
+       // evaluate mean multiplicity in transverse regions
+       Double_t err1=0.;
+       Double_t err2=0.;
+       Double_t err3=0.;
+       Double_t mean=(tmpetaphi[i]->IntegralAndError(1,5,err1)+tmpetaphi[i]->IntegralAndError(27,36,err2)+tmpetaphi[i]->IntegralAndError(58,62,err3))/(20.);
+       err1=TMath::Sqrt(err1*err1+err2*err2+err3*err3)/20.;
+       Printf("Branch: %d  MeanTransvMult: %f err: %f",i,mean,err1);
+        }
+  legd->Draw("same");
+
+  //8) TRACK-pT DISTRIBUTION IN FULL REGION
+  TCanvas *c7=((TCanvas*)arrC->At(7));
+  c7->cd();
+  SetStyle();
+  gPad->SetLogy();
+  gPad->SetLogx();
+  gPad->SetTitle("Full region (projection X)");
+  TH1D *tmpfull[20];
+  for (Int_t i= 0; i<br; i++){
+       tmpfull[i]=new TH1D();
+       tmpfull[i]=hptfull[i]->ProjectionX(Form("full%d",i),binstartprojection);
+       tmpfull[i]->GetYaxis()->SetTitle("1/n_{lj} dN/dp_{T}");
+        tmpfull[i]->Scale(1./(tmpfull[i]->GetBinWidth(1)*normphi[i]));
+       tmpfull[i]->SetMarkerStyle(20);
+       tmpfull[i]->SetMarkerColor(i+1);
+       tmpfull[i]->SetMaximum(100.);
+       if (i==0)tmpfull[i]->Draw("E1");
+       else tmpfull[i]->Draw("same"); 
+        }
+  legp->Draw("same");
+
+  //9) TRACK-pT DISTRIBUTION IN TRANSVERSE (MIN+MAX) REGION
+  TCanvas *c8=((TCanvas*)arrC->At(8));
+  c8->cd();
+  SetStyle();
+  gPad->SetLogy();
+  gPad->SetLogx();
+  gPad->SetTitle("Transverse regions (projection X)");
+  TH1D *tmptransv[20];
+  for (Int_t i= 0; i<br; i++){
+       tmptransv[i]=new TH1D();
+       tmptransv[i]=hpttransv[i]->ProjectionX(Form("transv%d",i),binstartprojection);
+       tmptransv[i]->GetYaxis()->SetTitle("1/n_{lj} dN/dp_{T}");
+        tmptransv[i]->Scale(1./(tmptransv[i]->GetBinWidth(1)*normphi[i]));
+       tmptransv[i]->SetMarkerStyle(20);
+       tmptransv[i]->SetMarkerColor(i+1);
+       tmptransv[i]->SetMaximum(10.);
+       if (i==0)tmptransv[i]->Draw("E1");
+       else tmptransv[i]->Draw("same"); 
+        }
+  legp->Draw("same");
+}
diff --git a/PWG4/JetTasks/AliHistogramsUE.h b/PWG4/JetTasks/AliHistogramsUE.h
new file mode 100644 (file)
index 0000000..3fdc1c9
--- /dev/null
@@ -0,0 +1,102 @@
+#ifndef ALIHISTOGRAMSUE_H
+#define ALIHISTOGRAMSUE_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+////////////////////////////////////////////////
+//--------------------------------------------- 
+// Class  to handle histograms for UE analysis
+//---------------------------------------------
+////////////////////////////////////////////////
+
+class TH1F;
+class TH2F;
+class TH1I;
+class TObjArray;
+class TProfile;
+class TTree;
+class TVector3;
+
+class  AliHistogramsUE : public TObject
+  {
+  public:
+    AliHistogramsUE();
+    AliHistogramsUE(TList * list);
+    virtual           ~AliHistogramsUE() { }
+    AliHistogramsUE(const  AliHistogramsUE &det);
+    AliHistogramsUE&   operator=(const  AliHistogramsUE &det);
+    
+    TObjArray*     CreateCanvas(const Int_t ncanv);
+    TObjArray*     GetHistosForPlotting(TString file, TString branches);
+    TList*        CreateHistos(Int_t bins, Double_t min, Double_t max, Double_t etacut);
+    void           DrawUE(Int_t debug);  //to draw final plots (normalized)
+    void           FillHistogram(const char* name,Double_t fillX); //One dimensional
+    void           FillHistogram(const char* name,Int_t fillX); //One dimensional
+    void           FillHistogram(const char* name,Double_t fillX, Double_t fillY); //Two dimensional
+    void           FillHistogram(const char* name,Double_t fillX, Double_t fillY, Double_t weight); //Two dimensional
+    void           FillHistogram(const char* name,Double_t fillX, Int_t fillY, Double_t weight); //Two dimensional
+    TList*        GetListOfHistos();
+    TH1F*            GetTrials()       {return fh1Trials;}
+    TProfile*        GetXsec()         {return fh1Xsec;}
+    void           PlotBranchesUE(TString file, TString branches, Double_t minJetProjection);  //TO BE CALLED BY EXTERNAL MACRO !!!             
+    void           SetStyle(); 
+  protected:
+
+  private:
+    
+    Int_t          fBinsPtInHist;    // Number of pT bins in histograms
+    Double_t       fMinJetPtInHist;  // Minimum jet pT in histograms
+    Double_t       fMaxJetPtInHist;  // Maximum jet pT in histograms
+    Double_t       fTrackEtaCut;     // Track eta cut  
+    TList*         fListOfHistos;    //  Output list of histograms
+    
+   
+    // Histograms
+    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; //!
+    
+    TH2F*  fhRegForwardMult;         //!
+    TH2F*  fhRegForwardSumPtvsMult;  //!
+    TH2F*  fhRegBackwardMult;        //!
+    TH2F*  fhRegBackwardSumPtvsMult; //!
+    TH2F*  fhRegForwardPartPtDistVsEt; //!
+    TH2F*  fhRegBackwardPartPtDistVsEt; //!
+    TH2F*  fhRegTransMult;         //!
+    TH2F*  fhRegTransSumPtVsMult;    //!
+    TH2F*  fhMinRegSumPtJetPtBin;    //!
+    TH2F*  fhMaxRegSumPtJetPtBin;    //!
+    TH1F*  fhVertexMult;             //!
+    TProfile*  fh1Xsec;                    //!         
+    TH1F*  fh1Trials;               //!
+
+    ClassDef( AliHistogramsUE, 1 ); // Class to manage histograms in UE analysis
+  };
+
+#endif
+
+    
diff --git a/PWG4/JetTasks/AliPWG4CosmicCandidates.cxx b/PWG4/JetTasks/AliPWG4CosmicCandidates.cxx
new file mode 100644 (file)
index 0000000..9b374bf
--- /dev/null
@@ -0,0 +1,356 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------------
+// Analysis task looking for cosmic candidates
+// Task checks if particles are back-to-back in eta and phi
+//             
+//
+// Author : Marta Verweij - UU - marta.verweij@cern.ch
+//-----------------------------------------------------------------------
+
+#include "TVector3.h"
+#include <iostream>
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TList.h"
+#include "TChain.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+
+#include "AliESDtrack.h"
+#include "AliESDtrackCuts.h"
+#include "AliExternalTrackParam.h"
+
+#include "AliLog.h"
+
+#include "AliPWG4CosmicCandidates.h"
+
+
+using namespace std;
+ClassImp(AliPWG4CosmicCandidates)
+
+//using namespace std; //required for resolving the 'cout' symbol
+
+//________________________________________________________________________
+AliPWG4CosmicCandidates::AliPWG4CosmicCandidates()
+  : AliAnalysisTaskSE(),
+    fTrackCuts(0), 
+    fPtMin(5.),
+    fMaxCosmicAngle(0.002),
+    fPtSignedCosmicCandidates(0),
+    fDeltaPtCosmicCandidates(0),
+    fDeltaPhiSumEta(0),
+    fDCAZCosmicCandidates(0),
+    fDCARCosmicCandidates(0),
+    fTheta(0),
+    fThetaZoom(0),
+    fThetaPt1Pt2(0),
+    fThetaPt1Pt2Signed(0),
+    fDeltaPhiSumEtaPt1(0),
+    fDeltaPhiSumEtaPt2(0),
+    fThetaDCAZ1DCAZ2(0),
+    fRisol(0),
+    fRisolTheta(0),
+    fHistListCosmics(0)
+{
+  // Default constructor
+}
+
+//________________________________________________________________________
+AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const char *name)
+  : AliAnalysisTaskSE(name),
+    fTrackCuts(0), 
+    fPtMin(5.),
+    fMaxCosmicAngle(0.002),
+    fPtSignedCosmicCandidates(0),
+    fDeltaPtCosmicCandidates(0),
+    fDeltaPhiSumEta(0),
+    fDCAZCosmicCandidates(0),
+    fDCARCosmicCandidates(0),
+    fTheta(0),
+    fThetaZoom(0),
+    fThetaPt1Pt2(0),
+    fThetaPt1Pt2Signed(0),
+    fDeltaPhiSumEtaPt1(0),
+    fDeltaPhiSumEtaPt2(0),
+    fThetaDCAZ1DCAZ2(0),
+    fRisol(0),
+    fRisolTheta(0),
+    fHistListCosmics(0)
+{
+  // Constructor. Initialization of Inputs and Outputs
+
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  //  DefineInput(0, TChain::Class());
+  // Output slot #0 id reserved by the base class for AOD
+  // Output slot #1 writes into a TList
+  DefineOutput(1, TList::Class());
+  // Output slot #2 writes into a AliESDtrackCuts
+  DefineOutput(2, AliESDtrackCuts::Class());
+}
+//________________________________________________________________________
+AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const AliPWG4CosmicCandidates &res)
+  : AliAnalysisTaskSE(res),
+    fTrackCuts(0), 
+    fPtMin(5.),
+    fMaxCosmicAngle(0.002),
+    fPtSignedCosmicCandidates(0),
+    fDeltaPtCosmicCandidates(0),
+    fDeltaPhiSumEta(0),
+    fDCAZCosmicCandidates(0),
+    fDCARCosmicCandidates(0),
+    fTheta(0),
+    fThetaZoom(0),
+    fThetaPt1Pt2(0),
+    fThetaPt1Pt2Signed(0),
+    fDeltaPhiSumEtaPt1(0),
+    fDeltaPhiSumEtaPt2(0),
+    fThetaDCAZ1DCAZ2(0),
+    fRisol(0),
+    fRisolTheta(0),
+    fHistListCosmics(0)
+{
+    // Dummy copy constructor
+}
+
+//________________________________________________________________________
+AliPWG4CosmicCandidates& AliPWG4CosmicCandidates::operator=(const AliPWG4CosmicCandidates& /*trclass*/)
+{
+    // Dummy assignment operator
+    return *this;
+}
+
+//________________________________________________________________________
+void AliPWG4CosmicCandidates::LocalInit()
+{
+  //
+  // Only called once at beginning
+  //
+  PostData(2,fTrackCuts);
+}
+
+//________________________________________________________________________
+void AliPWG4CosmicCandidates::UserCreateOutputObjects()
+{
+  //Create output objects
+  // Called once
+  AliDebug(2,Form(">> AliPWG4CosmicCandidates::UserCreateOutputObjects \n")); 
+
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE); 
+  
+  OpenFile(1);
+  fHistListCosmics = new TList();
+
+  fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
+  fHistListCosmics->Add(fNEventAll);
+  fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
+  fHistListCosmics->Add(fNEventSel);
+
+  Float_t fgkPtMin=0.;
+  Float_t fgkPtMax=100.;
+  Int_t fgkNPtBins= (int)(fgkPtMax-fgkPtMin);
+
+  Int_t fgkNPhiBins=18;
+  Float_t kMinPhi = -0.5*TMath::Pi();
+  Float_t kMaxPhi = 3./2.*TMath::Pi();
+
+  Int_t fgkNThetaBins=fgkNPhiBins*8;
+  Float_t kMinTheta = -0.5*TMath::Pi();
+  Float_t kMaxTheta = 3./2.*TMath::Pi();
+
+  Int_t fgkNDCARBins=80;
+  Float_t fgkDCARMin = -0.2;
+  Float_t fgkDCARMax = 0.2;
+  Double_t *binsDCAR=new Double_t[fgkNDCARBins+1];
+  for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Double_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Double_t)i ;
+
+  Int_t fgkNDCAZBins=80;
+  Float_t fgkDCAZMin = -2.;
+  Float_t fgkDCAZMax = 2.;
+  Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
+  for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
+
+  fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*(int)(fgkPtMax-fgkPtMin), -1.*fgkPtMax, fgkPtMax);
+  fHistListCosmics->Add(fPtSignedCosmicCandidates);  
+
+  fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
+  fHistListCosmics->Add(fDeltaPtCosmicCandidates);  
+
+  fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
+  fHistListCosmics->Add(fDeltaPhiSumEta);  
+
+  fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
+  fHistListCosmics->Add(fDCAZCosmicCandidates);
+
+  fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
+  fHistListCosmics->Add(fDCARCosmicCandidates);
+
+  fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
+  fHistListCosmics->Add(fTheta);
+
+  fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
+  fHistListCosmics->Add(fThetaZoom);
+
+  fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
+  fHistListCosmics->Add(fThetaPt1Pt2);
+
+  fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,4*(int)(fgkPtMax-fgkPtMin),-1.*fgkPtMax,fgkPtMax,4*(int)(fgkPtMax-fgkPtMin),-1.*fgkPtMax,fgkPtMax);
+  fHistListCosmics->Add(fThetaPt1Pt2Signed);
+
+  fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
+  fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
+
+  fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
+  fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
+
+  fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
+  fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
+
+  fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
+  fHistListCosmics->Add(fRisol);
+
+  fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
+  fHistListCosmics->Add(fRisolTheta);
+
+  TH1::AddDirectory(oldStatus);   
+
+  if(binsDCAR) delete [] binsDCAR;
+  if(binsDCAZ) delete [] binsDCAZ;
+
+}
+
+//________________________________________________________________________
+void AliPWG4CosmicCandidates::UserExec(Option_t *) 
+{
+//   // Main loop
+//   // Called for each event
+
+  // All events without selection
+  fNEventAll->Fill(0.);
+   
+  if (!fInputEvent) {
+    AliDebug(2,Form("ERROR: fESD not available"));
+    cout << "ERROR: fESD not available" << endl;
+    PostData(1, fHistListCosmics);
+    return;
+  }
+
+  const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
+  // Need vertex cut
+  if (!vtx || vtx->GetNContributors() < 2) {
+    // Post output data
+    PostData(1, fHistListCosmics);
+    return;
+  }
+
+  //  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
+  double primVtx[3];
+  vtx->GetXYZ(primVtx);
+  if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
+    // Post output data
+    PostData(1, fHistListCosmics);
+    return;
+  }
+  if(!fInputEvent->GetNumberOfTracks() || fInputEvent->GetNumberOfTracks()<2){ 
+    // Post output data
+    PostData(1, fHistListCosmics);
+    return;
+  }
+  Int_t nTracks = fInputEvent->GetNumberOfTracks();
+
+  if(!fTrackCuts) {
+   // Post output data
+    PostData(1, fHistListCosmics);
+    return;
+  }
+
+ fNEventSel->Fill(0.);
+
+  Float_t dcaR[2] = {0.,0.};
+  Float_t dcaZ[2] = {0.,0.};
+
+  // Track loop to fill a pT spectrum
+  for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
+
+    AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
+    if (!track1)  continue;
+    if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
+    if(track1->Pt()<fPtMin) continue;
+    //Start 2nd track loop to look for correlations
+    for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
+      AliESDtrack *track2 = (AliESDtrack*)fInputEvent->GetTrack(iTrack2);
+      if(!track2) continue;
+      if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
+         
+      //Check if back-to-back
+      Double_t mom1[3],mom2[3];
+      track1->GetPxPyPz(mom1);
+      track2->GetPxPyPz(mom2);
+      //     Double_t cosTheta = (mom1[0]*mom2[0]+mom1[1]*mom2[1]+mom1[2]*mom2[2])/( TMath::Sqrt(mom1[0]*mom1[0]+mom1[1]*mom1[1]+mom1[2]*mom1[2])*TMath::Sqrt(mom2[0]*mom2[0]+mom2[1]*mom2[1]+mom2[2]*mom2[2]) );
+      TVector3 momv1(mom1[0],mom1[1],mom1[2]);
+      TVector3 momv2(mom2[0],mom2[1],mom2[2]);
+      //Double_t theta = momv1.Angle(momv2);
+      Double_t theta = momv1.Phi()-momv2.Phi();
+      if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
+    
+      fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
+      Float_t deltaPhi = track1->Phi()-track2->Phi();
+      if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
+      fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
+
+      track1->GetImpactParameters(dcaR[0],dcaZ[0]);
+      track2->GetImpactParameters(dcaR[1],dcaZ[1]);
+
+      if(track2->Pt()<0.5) continue;
+      Double_t rIsol = TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) );
+      fRisol->Fill(rIsol); //Fill R histogram
+      if(track2->Pt()<fPtMin) continue;
+
+       fTheta->Fill(theta);
+       fThetaZoom->Fill(theta);
+       fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
+       fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()*track1->Pt(),track2->GetSign()*track2->Pt());
+       fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
+       fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
+       fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
+       fRisolTheta->Fill(rIsol,theta);
+       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
+         fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
+         fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
+         fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
+       }
+
+    } // track2 loop
+    
+  } // track1 loop 
+
+  PostData(1,fHistListCosmics);
+}      
+
+//________________________________________________________________________
+void AliPWG4CosmicCandidates::Terminate(Option_t *) 
+{
+  // Called once at the end of the query
+
+  
+}
diff --git a/PWG4/JetTasks/AliPWG4CosmicCandidates.h b/PWG4/JetTasks/AliPWG4CosmicCandidates.h
new file mode 100644 (file)
index 0000000..0daada5
--- /dev/null
@@ -0,0 +1,69 @@
+#ifndef AliPWG4CosmicCandidates_cxx\r
+#define AliPWG4CosmicCandidates_cxx\r
+\r
+// Analysis task looking for cosmic candidates\r
+// Authors: Marta Verweij marta.verweij@cern.ch\r
+\r
+class TH1F;\r
+class TH2F;\r
+class TH3F;\r
+class TList;\r
+class AliESDEvent;\r
+class AliESDfriend;\r
+class AliESDfriendTrack;\r
+class AliMCEvent;\r
+class AliVEvent;\r
+class AliESDtrackCuts;\r
+class AliESDtrack;\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliPWG4CosmicCandidates : public AliAnalysisTaskSE {\r
+ public:\r
+  AliPWG4CosmicCandidates();\r
+  AliPWG4CosmicCandidates(const char *name);\r
+  AliPWG4CosmicCandidates(const AliPWG4CosmicCandidates &res);\r
+  AliPWG4CosmicCandidates& operator=(const AliPWG4CosmicCandidates& trclass);\r
+  virtual ~AliPWG4CosmicCandidates() {;}\r
+\r
+  virtual void   LocalInit();\r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   Terminate(Option_t *);\r
+  \r
+  //Setters\r
+  void SetCuts(AliESDtrackCuts* trackCuts) {fTrackCuts = trackCuts;}\r
+  void SetPtMin(Double_t ptmin)            {fPtMin = ptmin;}\r
+  void SetMaxCosmicAngle(Double_t angle)   {fMaxCosmicAngle = angle;}\r
+\r
+ private:\r
+  AliESDtrackCuts *fTrackCuts;    // Standard trackCuts for global tracks\r
+\r
+  Double_t fPtMin;                // Minimal pt for cosmic candidate \r
+  Double_t fMaxCosmicAngle;       // Max deviation from pi (angle between two tracks) in case of cosmic candidate\r
+\r
+  TH1F *fNEventAll;                             //! Event counter\r
+  TH1F *fNEventSel;                             //! Event counter: Selected events for analysis\r
+  \r
+  TH1F *fPtSignedCosmicCandidates;              //! Cosmic Candidates\r
+  TH1F *fDeltaPtCosmicCandidates;               //! Cosmic Candidates Delta Pt\r
+  TH2F *fDeltaPhiSumEta;                        //! Cosmic Candidates Delta Phi vs Sum Eta\r
+  TH2F *fDCAZCosmicCandidates;                  //! Cosmic Candidates DCAZ track1 vs track2\r
+  TH2F *fDCARCosmicCandidates;                  //! Cosmic Candidates DCAR track1 vs track2\r
+  TH1F *fTheta;                                 //! Angle \theta between cosmic candidates in 3D space\r
+  TH1F *fThetaZoom;                             //! Angle between cosmic candidates in 3D space zoomed into back-to-back region\r
+  TH3F *fThetaPt1Pt2;                           //! Angle theta vs Pt1 vs Pt2\r
+  TH3F *fThetaPt1Pt2Signed;                     //! Angle theta vs Pt1 vs Pt2\r
+  TH3F *fDeltaPhiSumEtaPt1;                     //! Delta Phi vs Sum Eta vs Pt1\r
+  TH3F *fDeltaPhiSumEtaPt2;                     //! Delta Phi vs Sum Eta vs Pt2\r
+  TH3F *fThetaDCAZ1DCAZ2;                       //! Angle theta vs DCAZ1 vs DCAZ2\r
+  TH1F *fRisol;                                 //! Isolation R\r
+  TH2F *fRisolTheta;                            //! Isolation R vs Theta\r
+\r
+  TList *fHistListCosmics;                      //! List of Histograms for cosmic candidates  \r
+\r
+  \r
+  ClassDef(AliPWG4CosmicCandidates, 1);\r
+};\r
+\r
+#endif\r
index 870f46f..0ebda15 100644 (file)
@@ -74,7 +74,6 @@ AliPWG4HighPtQATPConly::AliPWG4HighPtQATPConly(): AliAnalysisTask("AliPWG4HighPt
   fPtAllminPtTPCvsPtAllEtaPos(0),
   fPtAllminPtTPCvsPtAllEtaNeg(0),
   fPtAllminPtTPCvsPtAllNPointTPC(0),
-  fPtAllminPtTPCvsPtAllNPointTPCS(0),
   fPtAllminPtTPCvsPtAllDCAR(0),
   fPtAllminPtTPCvsPtAllDCAZ(0),
   fPtAllminPtTPCvsPtAllPhi(0),
@@ -93,7 +92,6 @@ AliPWG4HighPtQATPConly::AliPWG4HighPtQATPConly(): AliAnalysisTask("AliPWG4HighPt
   fPtITSouterminPtTPCvsPtAllEtaPos(0),
   fPtITSouterminPtTPCvsPtAllEtaNeg(0),
   fPtITSouterminPtTPCvsPtAllNPointTPC(0),
-  fPtITSouterminPtTPCvsPtAllNPointTPCS(0),
   fPtITSouterminPtTPCvsPtAllDCAR(0),
   fPtITSouterminPtTPCvsPtAllDCAZ(0),
   fPtITSouterminPtTPCvsPtAllPhi(0),
@@ -131,7 +129,6 @@ AliPWG4HighPtQATPConly::AliPWG4HighPtQATPConly(): AliAnalysisTask("AliPWG4HighPt
   fPtITSminPtTPCvsPtITSEtaPos(0),
   fPtITSminPtTPCvsPtITSEtaNeg(0),
   fPtITSminPtTPCvsPtITSNPointTPC(0),
-  fPtITSminPtTPCvsPtITSNPointTPCS(0),
   fPtITSminPtTPCvsPtITSDCAR(0),
   fPtITSminPtTPCvsPtITSDCAZ(0),
   fPtITSminPtTPCvsPtITSPhi(0),
@@ -145,9 +142,6 @@ AliPWG4HighPtQATPConly::AliPWG4HighPtQATPConly(): AliAnalysisTask("AliPWG4HighPt
   fPtITSminPtTPCvsNPointITSPhi(0),
   fPtITSminPtTPCvsRel1PtUncertaintyPhi(0),
   fPtRel1PtUncertaintyChi2PerClusTPC(0),
-  fPtNPointTPCSChi2PerClusTPC(0),
-  fPtNPointTPCSRel1PtUncertainty(0),
-  fPtRel1PtUncertaintyChi2PerClusTPCSharedSel(0),
   fHistListITS(0),
   fPtSignedCosmicCandidates(0),
   fDeltaPtCosmicCandidates(0),
@@ -183,7 +177,6 @@ AliPWG4HighPtQATPConly::AliPWG4HighPtQATPConly(const char *name):
   fPtAllminPtTPCvsPtAllEtaPos(0),
   fPtAllminPtTPCvsPtAllEtaNeg(0),
   fPtAllminPtTPCvsPtAllNPointTPC(0),
-  fPtAllminPtTPCvsPtAllNPointTPCS(0),
   fPtAllminPtTPCvsPtAllDCAR(0),
   fPtAllminPtTPCvsPtAllDCAZ(0),
   fPtAllminPtTPCvsPtAllPhi(0),
@@ -202,7 +195,6 @@ AliPWG4HighPtQATPConly::AliPWG4HighPtQATPConly(const char *name):
   fPtITSouterminPtTPCvsPtAllEtaPos(0),
   fPtITSouterminPtTPCvsPtAllEtaNeg(0),
   fPtITSouterminPtTPCvsPtAllNPointTPC(0),
-  fPtITSouterminPtTPCvsPtAllNPointTPCS(0),
   fPtITSouterminPtTPCvsPtAllDCAR(0),
   fPtITSouterminPtTPCvsPtAllDCAZ(0),
   fPtITSouterminPtTPCvsPtAllPhi(0),
@@ -240,7 +232,6 @@ AliPWG4HighPtQATPConly::AliPWG4HighPtQATPConly(const char *name):
   fPtITSminPtTPCvsPtITSEtaPos(0),
   fPtITSminPtTPCvsPtITSEtaNeg(0),
   fPtITSminPtTPCvsPtITSNPointTPC(0),
-  fPtITSminPtTPCvsPtITSNPointTPCS(0),
   fPtITSminPtTPCvsPtITSDCAR(0),
   fPtITSminPtTPCvsPtITSDCAZ(0),
   fPtITSminPtTPCvsPtITSPhi(0),
@@ -254,9 +245,6 @@ AliPWG4HighPtQATPConly::AliPWG4HighPtQATPConly(const char *name):
   fPtITSminPtTPCvsNPointITSPhi(0),
   fPtITSminPtTPCvsRel1PtUncertaintyPhi(0),
   fPtRel1PtUncertaintyChi2PerClusTPC(0),
-  fPtNPointTPCSChi2PerClusTPC(0),
-  fPtNPointTPCSRel1PtUncertainty(0),
-  fPtRel1PtUncertaintyChi2PerClusTPCSharedSel(0),
   fHistListITS(0),
   fPtSignedCosmicCandidates(0),
   fDeltaPtCosmicCandidates(0),
@@ -415,12 +403,6 @@ void AliPWG4HighPtQATPConly::CreateOutputObjects() {
   Double_t *binsNPointTPC=new Double_t[fgkNNPointTPCBins+1];
   for(Int_t i=0; i<=fgkNNPointTPCBins; i++) binsNPointTPC[i]=(Double_t)fgkNPointTPCMin + (fgkNPointTPCMax-fgkNPointTPCMin)/fgkNNPointTPCBins*(Double_t)i ;
 
-  Int_t fgkNNPointTPCSBins=50;
-  Float_t fgkNPointTPCSMin = 0.;
-  Float_t fgkNPointTPCSMax = 1.;
-  Double_t *binsNPointTPCS=new Double_t[fgkNNPointTPCSBins+1];
-  for(Int_t i=0; i<=fgkNNPointTPCSBins; i++) binsNPointTPCS[i]=(Double_t)fgkNPointTPCSMin + (fgkNPointTPCSMax-fgkNPointTPCSMin)/fgkNNPointTPCSBins*(Double_t)i ;
-
   Int_t fgkNDCARBins=80;
   Float_t fgkDCARMin = -0.2;
   Float_t fgkDCARMax = 0.2;
@@ -489,12 +471,6 @@ void AliPWG4HighPtQATPConly::CreateOutputObjects() {
   fPtAllminPtTPCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
   fHistList->Add(fPtAllminPtTPCvsPtAllNPointTPC);
 
-  fPtAllminPtTPCvsPtAllNPointTPCS = new TH3F("fPtAllminPtTPCvsPtAllNPointTPCS","PtAllminPtTPCvsPtAllNPointTPCS",fgkNPtBins, binsPt,fgkNResPtBins,binsResPt,fgkNNPointTPCSBins,binsNPointTPCS);
-  fPtAllminPtTPCvsPtAllNPointTPCS->SetXTitle("p_{t}^{Global}");
-  fPtAllminPtTPCvsPtAllNPointTPCS->SetYTitle("(1/p_{t}^{Global}-1/p_{t}^{TPC})/(1/p_{t}^{Global})");
-  fPtAllminPtTPCvsPtAllNPointTPCS->SetZTitle("N_{point,TPC}^{Shared}/N_{point,TPC}");
-  fHistList->Add(fPtAllminPtTPCvsPtAllNPointTPCS);
-
   fPtAllminPtTPCvsPtAllDCAR = new TH3F("fPtAllminPtTPCvsPtAllDCAR","PtAllminPtTPCvsPtAllDCAR",fgkNPtBins, binsPt,fgkNResPtBins,binsResPt,fgkNDCARBins,binsDCAR);
   fPtAllminPtTPCvsPtAllDCAR->SetXTitle("p_{t}^{Global}");
   fPtAllminPtTPCvsPtAllDCAR->SetYTitle("(1/p_{t}^{Global}-1/p_{t}^{TPC})/(1/p_{t}^{Global})");
@@ -597,12 +573,6 @@ void AliPWG4HighPtQATPConly::CreateOutputObjects() {
   fPtITSouterminPtTPCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
   fHistList->Add(fPtITSouterminPtTPCvsPtAllNPointTPC);
 
-  fPtITSouterminPtTPCvsPtAllNPointTPCS = new TH3F("fPtITSouterminPtTPCvsPtAllNPointTPCS","PtITSouterminPtTPCvsPtAllNPointTPCS",fgkNPtBins, binsPt,fgkNResPtBins,binsResPt,fgkNNPointTPCSBins,binsNPointTPCS);
-  fPtITSouterminPtTPCvsPtAllNPointTPCS->SetXTitle("p_{t}^{Global}");
-  fPtITSouterminPtTPCvsPtAllNPointTPCS->SetYTitle("(1/p_{t}^{Global}-1/p_{t}^{TPC})/(1/p_{t}^{Global})");
-  fPtITSouterminPtTPCvsPtAllNPointTPCS->SetZTitle("N_{point,TPC}^{Shared}/N_{point,TPC}");
-  fHistList->Add(fPtITSouterminPtTPCvsPtAllNPointTPCS);
-
   fPtITSouterminPtTPCvsPtAllDCAR = new TH3F("fPtITSouterminPtTPCvsPtAllDCAR","PtITSouterminPtTPCvsPtAllDCAR",fgkNPtBins, binsPt,fgkNResPtBins,binsResPt,fgkNDCARBins,binsDCAR);
   fPtITSouterminPtTPCvsPtAllDCAR->SetXTitle("p_{t}^{Global}");
   fPtITSouterminPtTPCvsPtAllDCAR->SetYTitle("(1/p_{t}^{ITSouter}-1/p_{t}^{TPCinner})/(1/p_{t}^{ITSouter})");
@@ -784,12 +754,6 @@ void AliPWG4HighPtQATPConly::CreateOutputObjects() {
   fPtITSminPtTPCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
   fHistListITS->Add(fPtITSminPtTPCvsPtITSNPointTPC);
 
-  fPtITSminPtTPCvsPtITSNPointTPCS = new TH3F("fPtITSminPtTPCvsPtITSNPointTPCS","PtITSminPtTPCvsPtITSNPointTPCS",fgkNPtBins, binsPt,fgkNResPtBins,binsResPt,fgkNNPointTPCSBins,binsNPointTPCS);
-  fPtITSminPtTPCvsPtITSNPointTPCS->SetXTitle("p_{t}^{Global}");
-  fPtITSminPtTPCvsPtITSNPointTPCS->SetYTitle("(1/p_{t}^{Global}-1/p_{t}^{TPC})/(1/p_{t}^{Global})");
-  fPtITSminPtTPCvsPtITSNPointTPCS->SetZTitle("N_{point,TPC}^{Shared}/N_{point,TPC}");
-  fHistListITS->Add(fPtITSminPtTPCvsPtITSNPointTPCS);    
-
   fPtITSminPtTPCvsPtITSDCAR = new TH3F("fPtITSminPtTPCvsPtITSDCAR","PtITSminPtTPCvsPtITSDCAR",fgkNPtBins, binsPt,fgkNResPtBins,binsResPt,fgkNDCARBins,binsDCAR);
   fPtITSminPtTPCvsPtITSDCAR->SetXTitle("p_{t}^{Global}");
   fPtITSminPtTPCvsPtITSDCAR->SetYTitle("(1/p_{t}^{Global}-1/p_{t}^{TPC})/(1/p_{t}^{Global})");
@@ -868,24 +832,6 @@ void AliPWG4HighPtQATPConly::CreateOutputObjects() {
   fPtRel1PtUncertaintyChi2PerClusTPC->SetZTitle("#chi^{2}/(2*N_{clusters}^{TPC}-5)");
   fHistListITS->Add(fPtRel1PtUncertaintyChi2PerClusTPC);
 
-  fPtNPointTPCSChi2PerClusTPC = new TH3F("fPtNPointTPCSChi2PerClusTPC","PtITSminPtTPCvsPtITSNPointTPCS",fgkNPtBins, binsPt,fgkNNPointTPCSBins,binsNPointTPCS,fgkNChi2PerClusBins,binsChi2PerClus);
-  fPtNPointTPCSChi2PerClusTPC->SetXTitle("p_{t}^{global}");
-  fPtNPointTPCSChi2PerClusTPC->SetYTitle("N_{Point,TPC}^{Shared}/N_{Point,TPC}");
-  fPtNPointTPCSChi2PerClusTPC->SetZTitle("#chi^{2}/(2*N_{clusters}^{TPC}-5)");
-  fHistListITS->Add(fPtNPointTPCSChi2PerClusTPC);
-
-  fPtNPointTPCSRel1PtUncertainty = new TH3F("fPtNPointTPCSRel1PtUncertainty","PtITSminPtTPCvsPtITSNPointTPCS",fgkNPtBins, binsPt,fgkNNPointTPCSBins,binsNPointTPCS,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
-  fPtNPointTPCSRel1PtUncertainty->SetXTitle("p_{t}^{global}");
-  fPtNPointTPCSRel1PtUncertainty->SetYTitle("N_{Point,TPC}^{Shared}/N_{Point,TPC}");
-  fPtNPointTPCSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
-  fHistListITS->Add(fPtNPointTPCSRel1PtUncertainty);
-
-  fPtRel1PtUncertaintyChi2PerClusTPCSharedSel = new TH3F("fPtRel1PtUncertaintyChi2PerClusTPCSharedSel","PtITSminPtTPCvsPtITSRel1PtUncertainty",fgkNPtBins, binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty,fgkNChi2PerClusBins,binsChi2PerClus);
-  fPtRel1PtUncertaintyChi2PerClusTPCSharedSel->SetXTitle("p_{t}^{global}");
-  fPtRel1PtUncertaintyChi2PerClusTPCSharedSel->SetYTitle("Rel1PtUncertainty");
-  fPtRel1PtUncertaintyChi2PerClusTPCSharedSel->SetZTitle("#chi^{2}/(2*N_{clusters}^{TPC}-5)");
-  fHistListITS->Add(fPtRel1PtUncertaintyChi2PerClusTPCSharedSel);
-  
   fPtAllTPC = new TH1F("fPtAllTPC","PtAll",fgkNPtBins,binsPt);
   fHistListTPC->Add(fPtAllTPC);
   fPtSelTPC = new TH1F("fPtSelTPC","PtSel",fgkNPtBins,binsPt);
@@ -896,27 +842,31 @@ void AliPWG4HighPtQATPConly::CreateOutputObjects() {
   //****************************************************************************************************************//
   //                                              Cosmic Candidates                                                 //
   //****************************************************************************************************************//
+  Int_t fgkNThetaBins=fgkNPhiBins*8;
+  Float_t kMinTheta = -0.5*TMath::Pi();
+  Float_t kMaxTheta = 3./2.*TMath::Pi();
+
   fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*(int)(fgkPtMax-fgkPtMin), -1.*fgkPtMax, fgkPtMax);
   fHistListCosmics->Add(fPtSignedCosmicCandidates);  
   fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
   fHistListCosmics->Add(fDeltaPtCosmicCandidates);  
-  fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,0.,kMaxPhi,80, -2.,2.);
+  fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.);
   fHistListCosmics->Add(fDeltaPhiSumEta);  
   fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
   fHistListCosmics->Add(fDCAZCosmicCandidates);
   fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
   fHistListCosmics->Add(fDCARCosmicCandidates);
-  fTheta = new TH1F("fTheta","fTheta",fgkNPhiBins*8,0.,kMaxPhi);
+  fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
   fHistListCosmics->Add(fTheta);
   fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
   fHistListCosmics->Add(fThetaZoom);
-  fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNPhiBins*8,0.,kMaxPhi,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
+  fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
   fHistListCosmics->Add(fThetaPt1Pt2);
-  fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNPhiBins*4,0.,kMaxPhi,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
+  fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
   fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
-  fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNPhiBins*4,0.,kMaxPhi,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
+  fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
   fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
-  fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNPhiBins*8,0.,kMaxPhi,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
+  fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
   fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
 
   TH1::AddDirectory(oldStatus);   
@@ -926,7 +876,6 @@ void AliPWG4HighPtQATPConly::CreateOutputObjects() {
   if(binsPt) delete [] binsPt;
   if(binsResPt) delete [] binsResPt;
   if(binsNPointTPC) delete [] binsNPointTPC;
-  if(binsNPointTPCS) delete [] binsNPointTPCS;
   if(binsDCAR) delete [] binsDCAR;
   if(binsDCAZ) delete [] binsDCAZ;
   if(binsNPointITS) delete [] binsNPointITS;
@@ -1003,20 +952,13 @@ void AliPWG4HighPtQATPConly::Exec(Option_t *) {
   }
   
 
- const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
- if(vertex->GetNContributors()<1) {
-   // SPD vertex
-   vertex = fESD->GetPrimaryVertexSPD();
-   if(vertex->GetNContributors()<1) vertex = 0x0;
- }
-
   const AliESDVertex *vtx = fESD->GetPrimaryVertexTracks();
   // Need vertex cut
   if (vtx->GetNContributors() < 2) {
     // SPD vertex
     vtx = fESD->GetPrimaryVertexSPD();
     if(vtx->GetNContributors()<2) {
-      vertex = 0x0;
+      vtx = 0x0;
       // Post output data
       PostData(0, fHistList);
       PostData(1, fHistListTPC);
@@ -1116,7 +1058,6 @@ void AliPWG4HighPtQATPConly::Exec(Option_t *) {
       if(track->GetSign()<0.) fPtAllminPtTPCvsPtAllEtaNeg->Fill(pt,(1./pt-1./ptTPC)/(1./pt),track->Eta());
 
       fPtAllminPtTPCvsPtAllNPointTPC->Fill(pt,(1./pt-1./ptTPC)/(1./pt),nClustersTPC);
-      if(nClustersTPC>0.) fPtAllminPtTPCvsPtAllNPointTPCS->Fill(pt,(1./pt-1./ptTPC)/(1./pt),track->GetTPCnclsS()/nClustersTPC);
       fPtAllminPtTPCvsPtAllDCAR->Fill(pt,(1./pt-1./ptTPC)/(1./pt),dca2D);
       fPtAllminPtTPCvsPtAllDCAZ->Fill(pt,(1./pt-1./ptTPC)/(1./pt),dcaZ);
       fPtAllminPtTPCvsPtAllPhi->Fill(pt,(1./pt-1./ptTPC)/(1./pt),phi);
@@ -1145,7 +1086,6 @@ void AliPWG4HighPtQATPConly::Exec(Option_t *) {
          if(trackITSouter.GetSign()<0.) fPtITSouterminPtTPCvsPtAllEtaNeg->Fill(pt,(1./ptITSouter-1./ptTPC)/(1./ptITSouter),trackITSouter.Eta());
 
          fPtITSouterminPtTPCvsPtAllNPointTPC->Fill(pt,(1./ptITSouter-1./ptTPC)/(1./ptITSouter),nClustersTPC);
-         if(nClustersTPC>0.) fPtITSouterminPtTPCvsPtAllNPointTPCS->Fill(pt,(1./ptITSouter-1./ptTPC)/(1./ptITSouter),track->GetTPCnclsS()/nClustersTPC);
          fPtITSouterminPtTPCvsPtAllDCAR->Fill(pt,(1./ptITSouter-1./ptTPC)/(1./ptITSouter),dca2D);
          fPtITSouterminPtTPCvsPtAllDCAZ->Fill(pt,(1./ptITSouter-1./ptTPC)/(1./ptITSouter),dcaZ);
          fPtITSouterminPtTPCvsPtAllPhi->Fill(pt,(1./ptITSouter-1./ptTPC)/(1./ptITSouter),phi);
@@ -1206,7 +1146,6 @@ void AliPWG4HighPtQATPConly::Exec(Option_t *) {
       if(track->GetSign()>0.) fPtITSminPtTPCvsPtITSEtaPos->Fill(pt,(1./pt-1./ptTPC)/(1./pt),track->Eta());
       if(track->GetSign()<0.) fPtITSminPtTPCvsPtITSEtaNeg->Fill(pt,(1./pt-1./ptTPC)/(1./pt),track->Eta());
       fPtITSminPtTPCvsPtITSNPointTPC->Fill(pt,(1./pt-1./ptTPC)/(1./pt),nClustersTPC);
-      if(nClustersTPC>0.) fPtITSminPtTPCvsPtITSNPointTPCS->Fill(pt,(1./pt-1./ptTPC)/(1./pt),track->GetTPCnclsS()/nClustersTPC);
       fPtITSminPtTPCvsPtITSDCAR->Fill(pt,(1./pt-1./ptTPC)/(1./pt),dca2D);
       fPtITSminPtTPCvsPtITSDCAZ->Fill(pt,(1./pt-1./ptTPC)/(1./pt),dcaZ);
       fPtITSminPtTPCvsPtITSPhi->Fill(pt,(1./pt-1./ptTPC)/(1./pt),phi);
@@ -1221,9 +1160,6 @@ void AliPWG4HighPtQATPConly::Exec(Option_t *) {
       fPtITSminPtTPCvsRel1PtUncertaintyPhi->Fill((1./pt-1./ptTPC)/(1./pt),relUncertainty1Pt,phi);
 
       fPtRel1PtUncertaintyChi2PerClusTPC->Fill(pt,relUncertainty1Pt,chi2PerClusterTPC);
-      fPtNPointTPCSChi2PerClusTPC->Fill(pt,track->GetTPCnclsS()/nClustersTPC,chi2PerClusterTPC);
-      fPtNPointTPCSRel1PtUncertainty->Fill(pt,track->GetTPCnclsS()/nClustersTPC,relUncertainty1Pt);
-      if(track->GetTPCnclsS()/nClustersTPC>0.05) fPtRel1PtUncertaintyChi2PerClusTPCSharedSel->Fill(pt,relUncertainty1Pt,chi2PerClusterTPC);
     }//fTrackCutsITS loop
       
   }//ESD track loop
@@ -1262,18 +1198,20 @@ Bool_t AliPWG4HighPtQATPConly::IsCosmic(const AliESDtrack *track1 , Int_t trackN
     Double_t mom1[3],mom2[3];
     track1->GetPxPyPz(mom1);
     track2->GetPxPyPz(mom2);
-    Double_t cosTheta = (mom1[0]*mom2[0]+mom1[1]*mom2[1]+mom1[2]*mom2[2])/( TMath::Sqrt(mom1[0]*mom1[0]+mom1[1]*mom1[1]+mom1[2]*mom1[2])*TMath::Sqrt(mom2[0]*mom2[0]+mom2[1]*mom2[1]+mom2[2]*mom2[2]) );
-    Double_t theta = TMath::ACos(cosTheta);
-    if(cosTheta<0.) theta+=TMath::Pi();
-   //if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) { candidate1 = kTRUE; candidate2 = kTRUE;}
+    TVector3 momv1(mom1[0],mom1[1],mom1[2]);
+    TVector3 momv2(mom2[0],mom2[1],mom2[2]);
+    Double_t theta = momv1.Phi()-momv2.Phi();
+    if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
+
+    //if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) { candidate1 = kTRUE; candidate2 = kTRUE;}
     
-//    Double_t cosMaxCosmicAngle[2] = {TMath::Cos(TMath::Pi()-fMaxCosmicAngle),TMath::Cos(TMath::Pi()+fMaxCosmicAngle)};
-//    if(cosTheta >= cosMaxCosmicAngle[0] && cosTheta <= cosMaxCosmicAngle[1]) { 
+    //    Double_t cosMaxCosmicAngle[2] = {TMath::Cos(TMath::Pi()-fMaxCosmicAngle),TMath::Cos(TMath::Pi()+fMaxCosmicAngle)};
+    //    if(cosTheta >= cosMaxCosmicAngle[0] && cosTheta <= cosMaxCosmicAngle[1]) { 
     candidate1 = kTRUE; candidate2 = kTRUE;//}
     if(candidate2) {
       fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
       Float_t deltaPhi = track1->Phi()-track2->Phi();
-      if(deltaPhi<0.) deltaPhi+=2.*TMath::Pi();
+      if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
       fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
 
       track1->GetImpactParameters(dcaR[0],dcaZ[0]);
index a35d1f7..8b70748 100644 (file)
@@ -81,7 +81,6 @@ class AliPWG4HighPtQATPConly: public AliAnalysisTask {
   TH3F *fPtAllminPtTPCvsPtAllEtaPos;            //! Momentum resolution (global vs TPConly) vs Eta for positive particles
   TH3F *fPtAllminPtTPCvsPtAllEtaNeg;            //! Momentum resolution (global vs TPConly) vs Eta for negative particles
   TH3F *fPtAllminPtTPCvsPtAllNPointTPC;         //! Momentum resolution vs NPointTPC
-  TH3F *fPtAllminPtTPCvsPtAllNPointTPCS;        //! Momentum resolution vs NPointTPCShared/NPointTPC
   TH3F *fPtAllminPtTPCvsPtAllDCAR;              //! Momentum resolution vs DCAR
   TH3F *fPtAllminPtTPCvsPtAllDCAZ;              //! Momentum resolution vs DCAZ
   TH3F *fPtAllminPtTPCvsPtAllPhi;               //! Momentum resolution vs Phi
@@ -103,7 +102,6 @@ class AliPWG4HighPtQATPConly: public AliAnalysisTask {
   TH3F *fPtITSouterminPtTPCvsPtAllEtaPos;            //! Momentum resolution (global vs ITSouter-TPCinner) vs Eta for positive particles
   TH3F *fPtITSouterminPtTPCvsPtAllEtaNeg;            //! Momentum resolution (global vs ITSouter-TPCinner) vs Eta for negative particles
   TH3F *fPtITSouterminPtTPCvsPtAllNPointTPC;         //! Momentum resolution vs NPointTPC
-  TH3F *fPtITSouterminPtTPCvsPtAllNPointTPCS;        //! Momentum resolution vs NPointTPCS
   TH3F *fPtITSouterminPtTPCvsPtAllDCAR;              //! Momentum resolution vs DCAR
   TH3F *fPtITSouterminPtTPCvsPtAllDCAZ;              //! Momentum resolution vs DCAZ
   TH3F *fPtITSouterminPtTPCvsPtAllPhi;               //! Momentum resolution vs Phi
@@ -148,7 +146,6 @@ class AliPWG4HighPtQATPConly: public AliAnalysisTask {
   TH3F *fPtITSminPtTPCvsPtITSEtaPos;            //! Momentum resolution (global with ITSrefit vs TPConly) vs Eta for positive particles
   TH3F *fPtITSminPtTPCvsPtITSEtaNeg;            //! Momentum resolution (global with ITSrefit vs TPConly) vs Eta for negative particles
   TH3F *fPtITSminPtTPCvsPtITSNPointTPC;         //! Momentum resolution vs NPointTPC 
-  TH3F *fPtITSminPtTPCvsPtITSNPointTPCS;        //! Momentum resolution vs NPointTPCS 
   TH3F *fPtITSminPtTPCvsPtITSDCAR;              //! Momentum resolution vs DCAR
   TH3F *fPtITSminPtTPCvsPtITSDCAZ;              //! Momentum resolution vs DCAZ
   TH3F *fPtITSminPtTPCvsPtITSPhi;               //! Momentum resolution vs Phi
@@ -164,9 +161,6 @@ class AliPWG4HighPtQATPConly: public AliAnalysisTask {
   TH3F *fPtITSminPtTPCvsRel1PtUncertaintyPhi;   //! Momentum resolution vs Rel1PtUncertainty vs Phi
 
   TH3F *fPtRel1PtUncertaintyChi2PerClusTPC;     //! Global Pt vs relUncertainty1Pt vs Chi2PerClusTPC
-  TH3F *fPtNPointTPCSChi2PerClusTPC;            //! Global Pt vs NPointTPCShared/NPointTPC vs Chi2PerClusTPC
-  TH3F *fPtNPointTPCSRel1PtUncertainty;         //! Global Pt vs NPointTPCShared/NPointTPC vs relUncertainty1Pt
-  TH3F *fPtRel1PtUncertaintyChi2PerClusTPCSharedSel;     //! Global Pt vs relUncertainty1Pt vs Chi2PerClusTPC for NPointTPCShared/NPointTPC>0.05
 
   TList *fHistListITS; //! List of Histograms
 
index d73d3bf..50177e4 100644 (file)
@@ -65,7 +65,7 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpe
   fCFManagerNeg(0x0),
   fESD(0),
   fTrackCuts(0),
-  fTrigger(0),
+  fTrackCutsTPConly(0),
   fHistList(0),
   fNEventAll(0),
   fNEventSel(0)
@@ -82,7 +82,7 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
   fCFManagerNeg(0x0),
   fESD(0),
   fTrackCuts(),
-  fTrigger(0),
+  fTrackCutsTPConly(0),
   fHistList(0),
   fNEventAll(0),
   fNEventSel(0)
@@ -100,6 +100,7 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
   DefineOutput(2,AliCFContainer::Class());
   // Output slot #3 writes into a AliESDtrackCuts
   DefineOutput(3, AliESDtrackCuts::Class());
+  DefineOutput(4, AliESDtrackCuts::Class());
 }
 
 //________________________________________________________________________
@@ -109,55 +110,9 @@ void AliPWG4HighPtSpectra::LocalInit()
   // Only called once at beginning
   //
   PostData(3,fTrackCuts);
+  PostData(4,fTrackCutsTPConly);
 }
 
-// //___________________________________________________________________________
-// AliPWG4HighPtSpectra& AliPWG4HighPtSpectra::operator=(const AliPWG4HighPtSpectra& c) 
-// {
-//   //
-//   // Assignment operator
-//   //
-//   if (this!=&c) {
-//     AliAnalysisTask::operator=(c) ;
-//     fReadAODData = c.fReadAODData ;
-//     fCFManagerPos  = c.fCFManagerPos;
-//     fCFManagerNeg  = c.fCFManagerNeg;
-//     fHistList = c.fHistList;
-//     fNEventAll = c.fNEventAll;
-//     fNEventSel = c.fNEventSel;
-//   }
-//   return *this;
-// }
-
-// //___________________________________________________________________________
-// AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const AliPWG4HighPtSpectra& c) :
-//   AliAnalysisTask(c),
-//   fReadAODData(c.fReadAODData),
-//   fCFManagerPos(c.fCFManagerPos),
-//   fCFManagerNeg(c.fCFManagerNeg),
-//   fESD(c.fESD),
-//   fTrackCuts(c.fTrackCuts),
-//   fTrigger(c.fTrigger),
-//   fHistList(c.fHistList),
-//   fNEventAll(c.fNEventAll),
-//   fNEventSel(c.fNEventSel)
-// {
-//   //
-//   // Copy Constructor
-//   //
-// }
-
-// //___________________________________________________________________________
-// AliPWG4HighPtSpectra::~AliPWG4HighPtSpectra() {
-//   //
-//   //destructor
-//   //
-//   Info("~AliPWG4HighPtSpectra","Calling Destructor");
-//   if (fCFManagerPos)           delete fCFManagerPos ;
-//   if (fCFManagerNeg)           delete fCFManagerNeg ;
-//   if (fNEventAll)              delete fNEventAll ;
-//   if (fNEventSel)              delete fNEventSel ;
-// }
 //________________________________________________________________________
 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *) 
 {
@@ -307,24 +262,27 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
       containerInputRec[3] = dca2D;
       containerInputRec[4] = chi2PerClusterTPC;
 
+      //Store TPC Inner Params for TPConly tracks
       containerInputTPConly[0] = trackTPC->Pt();
       containerInputTPConly[1] = trackTPC->Phi();
       containerInputTPConly[2] = trackTPC->Eta();
       containerInputTPConly[3] = dca2DTPC/10.; //Divide by 10 in order to store in same container. Should be corrected back when looking at output.
       containerInputTPConly[4] = chi2PerClusterTPCIter1;//TPC;
 
-      if (fTrackCuts->AcceptTrack(track)) {
-       if(track->GetSign()>0.) {
-         fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
-         fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
-       }
-       if(track->GetSign()<0.) {
-         fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
-         fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
+      AliESDtrack* trackTPCESD = fTrackCutsTPConly->GetTPCOnlyTrack(fESD, iTrack);
+      if(trackTPCESD) {
+       if (fTrackCutsTPConly->AcceptTrack(trackTPCESD)) {
+         if(trackTPC->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
+         if(trackTPC->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
        }
-       
+      }
+
+      if (fTrackCuts->AcceptTrack(track)) {
+       if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
+       if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
+
        
-       //Only fill the secondary particle container if MC information is available
+       //Only fill the MC containers if MC information is available
        if(eventHandler) {
          Int_t label = TMath::Abs(track->GetLabel());
          TParticle *particle = stack->Particle(label) ;
@@ -336,13 +294,17 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
          containerInputMC[3] = 0.0;      
          containerInputMC[4] = 0.0;  
 
-         if(particle->GetPDG()->Charge()>0.) {
-           fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
-         }
-         if(particle->GetPDG()->Charge()<0.) {
-           fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
+         //Container with primaries
+         if(stack->IsPhysicalPrimary(label)) {
+           if(particle->GetPDG()->Charge()>0.) {
+             fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
+           }
+           if(particle->GetPDG()->Charge()<0.) {
+             fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
+           }
          }
 
+         //Container with secondaries
          if (!stack->IsPhysicalPrimary(label) ) {
            if(particle->GetPDG()->Charge()>0.) {
              fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
@@ -353,14 +315,14 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
          }
        }
        
-      }
+      }//trackCuts
       
-    }
+    }//track loop
   
 
   //Fill MC containters if particles are findable
   if(eventHandler) {
-    for(int iPart = 1; iPart<(mcEvent->GetNumberOfTracks()); iPart++)//stack->GetNprimary();
+    for(int iPart = 1; iPart<(mcEvent->GetNumberOfPrimaries()); iPart++)//stack->GetNprimary();
       {
        AliMCParticle *mcPart  = (AliMCParticle*)mcEvent->GetTrack(iPart);
        if(!mcPart) continue;
@@ -374,15 +336,17 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
 
        int counter;
        Float_t trackLengthTPC = 0.;
-       if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(3,mcPart)) {
-         fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance) ;
-         trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
-         if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
-       }
-       if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(3,mcPart)) {
-         fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance) ;
-         trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
-         if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
+       if(stack->IsPhysicalPrimary(iPart)) {
+         if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(3,mcPart)) {
+           fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance) ;
+           trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
+           if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
+         }
+         if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(3,mcPart)) {
+           fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance) ;
+           trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
+           if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
+         }
        }
       }
   }
index c8a431a..ebd0666 100644 (file)
@@ -64,8 +64,7 @@ class AliPWG4HighPtSpectra : public AliAnalysisTask {
   
   //AliESDtrackCuts setters
   void SetCuts(AliESDtrackCuts* trackCuts) {fTrackCuts = trackCuts;}
-  //Select the trigger
-  void SelectTrigger(Int_t trig) { fTrigger = trig; } 
+  void SetCutsTPConly(AliESDtrackCuts* trackCuts) {fTrackCutsTPConly = trackCuts;}
 
   // Data types
   Bool_t IsReadAODData()   const {return fReadAODData;}
@@ -76,21 +75,20 @@ class AliPWG4HighPtSpectra : public AliAnalysisTask {
   const AliCFManager  *fCFManagerPos    ;  // pointer to the CF manager for positive charged particles
   const AliCFManager  *fCFManagerNeg    ;  // pointer to the CF manager for negative charged particles
  
-  AliESDEvent *fESD;              //! ESD object
-  //AliESDtrackCuts options. Must be setted in AddTaskPWG4HighPtQAMC.C. They correspond with different steps in container.
-  AliESDtrackCuts *fTrackCuts;    // trackCuts applied
-  Int_t fTrigger;                 //Trigger flag as defined in AliAnalysisHelperJetTasks.h 
+  AliESDEvent *fESD;                     //! ESD object
+  //AliESDtrackCuts options. Must be setted in AddTaskPWG4HighPTSpectra.C. They correspond with different steps in container.
+  AliESDtrackCuts *fTrackCuts;           // trackCuts applied to global tracks
+  AliESDtrackCuts *fTrackCutsTPConly;    // trackCuts applied to TPConly tracks
 
  private:
- AliPWG4HighPtSpectra(const AliPWG4HighPtSpectra&);
- AliPWG4HighPtSpectra& operator=(const AliPWG4HighPtSpectra&);
-
+  AliPWG4HighPtSpectra(const AliPWG4HighPtSpectra&);
+  AliPWG4HighPtSpectra& operator=(const AliPWG4HighPtSpectra&);
 
   // Histograms
   //Number of events
   TList *fHistList;            //! List of output histograms
-  TH1F *fNEventAll;            //! Event counter
-  TH1F *fNEventSel;            //! Event counter: Selected events for analysis
+  TH1F  *fNEventAll;            //! Event counter
+  TH1F  *fNEventSel;            //! Event counter: Selected events for analysis
 
   ClassDef(AliPWG4HighPtSpectra,1);
 };
index a7ed63d..eda33c3 100644 (file)
@@ -5,6 +5,8 @@
 #pragma link off all functions;
 
 #pragma link C++ class AliAnalysisTaskUE+;
+#pragma link C++ class AliAnalyseUE+;
+#pragma link C++ class AliHistogramsUE+;
 #pragma link C++ class AliAnalysisTaskJetServices+;
 #pragma link C++ class AliAnalysisTaskJetSpectrum+;
 #pragma link C++ class AliAnalysisTaskJetSpectrum2+;
@@ -15,6 +17,7 @@
 #pragma link C++ class AliPWG4HighPtQATPConly+;
 #pragma link C++ class AliPWG4HighPtQAMC+;
 #pragma link C++ class AliPWG4HighPtSpectra+;
+#pragma link C++ class AliPWG4CosmicCandidates+;
 #pragma link C++ class AliAnalysisTaskPWG4PidDetEx+;
 #pragma link C++ class AliJetSpectrumUnfolding+;
 #pragma link C++ class AliAnalysisTaskJetChem+;
index 54267f1..88d14d9 100644 (file)
@@ -1,6 +1,6 @@
 #-*- Mode: Makefile -*-
 
-SRCS = JetTasks/AliAnalysisTaskUE.cxx  JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisTaskJetSpectrum2.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnalysisTaskJetServices.cxx JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx JetTasks/AliJetSpectrumUnfolding.cxx JetTasks/AliAnalysisTaskJFSystematics.cxx JetTasks/AliAnalysisTaskJetCorrections.cxx JetTasks/AliAnalysisTaskThreeJets.cxx JetTasks/AliPWG4HighPtQATPConly.cxx JetTasks/AliPWG4HighPtQAMC.cxx JetTasks/AliPWG4HighPtSpectra.cxx JetTasks/AliAnalysisTaskJetChem.cxx 
+SRCS = JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliHistogramsUE.cxx JetTasks/AliAnalyseUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisTaskJetSpectrum2.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnalysisTaskJetServices.cxx JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx JetTasks/AliJetSpectrumUnfolding.cxx JetTasks/AliAnalysisTaskJFSystematics.cxx JetTasks/AliAnalysisTaskJetCorrections.cxx JetTasks/AliAnalysisTaskThreeJets.cxx JetTasks/AliPWG4HighPtQATPConly.cxx JetTasks/AliPWG4HighPtQAMC.cxx JetTasks/AliPWG4HighPtSpectra.cxx JetTasks/AliPWG4CosmicCandidates.cxx JetTasks/AliAnalysisTaskJetChem.cxx 
 
 HDRS:= $(SRCS:.cxx=.h) 
 
diff --git a/PWG4/macros/AddTaskPWG4CosmicCandidates.C b/PWG4/macros/AddTaskPWG4CosmicCandidates.C
new file mode 100644 (file)
index 0000000..2c5f337
--- /dev/null
@@ -0,0 +1,67 @@
+//DEFINITION OF A FEW CONSTANTS
+const Double_t maxDeltaTheta = 0.01;
+const Double_t ptMin = 5.;
+
+AliPWG4CosmicCandidates* AddTaskPWG4CosmicCandidates(int cuts=1)
+{
+  // Creates HighPtQATPConly analysis task and adds it to the analysis manager.
+  
+  // A. Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("AddTaskPWG4CosmicCandidates", "No analysis manager to connect to.");
+    return NULL;
+  }  
+
+  // B. Check the analysis type using the event handlers connected to the analysis
+  //    manager. The availability of MC handler can also be checked here.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler()) {
+    ::Error("AddTaskPWG4CosmicCandidates", "This task requires an input event handler");
+    return NULL;
+  }  
+
+  TString inputDataType = mgr->GetInputEventHandler()->GetDataType(); 
+
+  // C. Create the task, add it to manager.
+  //===========================================================================
+  //CREATE THE  CUTS -----------------------------------------------
+  //Use AliESDtrackCuts
+  AliESDtrackCuts *trackCuts = new AliESDtrackCuts("AliESDtrackCuts","Standard Cuts");
+  if(cuts==1) trackCuts=trackCuts->GetStandardITSTPCTrackCuts2009(kTRUE);//Primary Track Selection
+  trackCuts->SetEtaRange(-0.9,0.9);
+  trackCuts->SetPtRange(0.15, 1e10);
+
+
+  //Create the task
+  cout << "Create the task AliPWG4CosmicCandidates" << endl;
+  AliPWG4CosmicCandidates *taskPWG4CC = new AliPWG4CosmicCandidates(Form("AliPWG4CosmicCandidates%d",cuts));
+  taskPWG4CC->SetCuts(trackCuts);
+  taskPWG4CC->SetMaxCosmicAngle(maxDeltaTheta);//0.008);
+  taskPWG4CC->SetPtMin(ptMin);
+  taskPWG4CC->SelectCollisionCandidates();
+
+  //Add task to manager
+  mgr->AddTask(taskPWG4CC);
+
+  // E. Create ONLY the output containers for the data produced by the task.
+  // Get and connect other common input/output containers via the manager as below
+  //==============================================================================
+  //  TString commonFileName = AliAnalysisManager::GetCommonFileName();
+  //  commonFileName += ":PWG4_CosmicCandidates"; 
+  AliAnalysisDataContainer *cout_hist1 = mgr->CreateContainer(Form("cosmic_hists%d",cuts), TList::Class(), AliAnalysisManager::kOutputContainer, Form("%s:PWG4_CosmicCandidates%d", AliAnalysisManager::GetCommonFileName(),cuts));
+  AliAnalysisDataContainer *cout_cuts0 = mgr->CreateContainer(Form("cosmic_cuts%d",cuts), AliESDtrackCuts::Class(), AliAnalysisManager::kParamContainer, Form("%s:PWG4_CosmicCandidates%d",AliAnalysisManager::GetCommonFileName(),cuts));
+  
+  
+  //Connect input containter to manager
+  mgr->ConnectInput(taskPWG4CC,0,mgr->GetCommonInputContainer());
+
+  //Connect output containers to manager
+  mgr->ConnectOutput(taskPWG4CC,1,cout_hist1);
+  mgr->ConnectOutput(taskPWG4CC,2,cout_cuts0);
+  
+  // Return task pointer at the end
+  return taskPWG4CC;
+}
index d9e44b7..5277025 100644 (file)
@@ -79,6 +79,7 @@ Int_t       iPWG4JetChem       = 0;      // Jet chemistry
 Int_t       iPWG4PtQAMC        = 0;      // Marta's QA tasks 
 Int_t       iPWG4PtSpectra     = 0;      // Marta's QA tasks 
 Int_t       iPWG4PtQATPC       = 0;      // Marta's QA tasks 
+Int_t       iPWG4PtCosmics     = 0;      // Marta's Cosmics Taks 
 Int_t       iPWG4ThreeJets     = 0;      // Sona's thrust task
 Int_t       iPWG4KMeans        = 0;      // Andreas' KMeans task 
 Int_t       iPWG4Cluster       = 0;      // CKB cluster task 
@@ -210,6 +211,7 @@ void AnalysisTrainPWG4Jets(const char *analysis_mode="local",
    printf(":: use PWG4 Pt QA MC       %d\n",iPWG4PtQAMC);
    printf(":: use PWG4 Pt Spectra     %d\n",iPWG4PtSpectra);
    printf(":: use PWG4 Pt QA TPC      %d\n",iPWG4PtQATPC);     
+   printf(":: use PWG4 Cosmics        %d\n",iPWG4Cosmics);     
    printf(":: use PWG4 Three Jets     %d\n",iPWG4ThreeJets);
    printf(":: use PWG4 KMeans         %d\n",iPWG4KMeans);
    printf(":: use PWG4 Cluster        %d\n",iPWG4Cluster);
@@ -468,6 +470,16 @@ void AnalysisTrainPWG4Jets(const char *analysis_mode="local",
  if (!taskQATPC) ::Warning("AnalysisTrainPWG4Jets", "AliAnalysisTaskQATPC cannot run for this train conditions - EXCLUDED");
    }
 
+   if(iPWG4Cosmics){
+     gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskPWG4CosmicCandidates.C");
+
+     AliPWG4CosmicCandidates *taskPWG4CosmicCandidates = AddTaskPWG4CosmicCandidates(0);
+     taskPWG4CosmicCandidates = AddTaskPWG4CosmicCandidates(1);
+
+     if (!taskPWG4CosmicCandidates) ::Warning("AnalysisTrainPWG4Jets", "AddTaskPWG4CosmicCandidates cannot run for this train conditions - EXCLUDED");
+   }
+
+
    if(iPWG4PtSpectra){
      gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskPWG4HighPtSpectra.C");
      AliPWG4HighPtSpectra *taskPtSpectra = AddTaskPWG4HighPtSpectra();
@@ -736,6 +748,9 @@ void CheckModuleFlags(const char *mode) {
       iPWG4PtQAMC        = 0;
       if( iPWG4PtQATPC)::Info("AnalysisTrainPWG4Jets.C::CheckModuleFlags", "PWG4 PtTPC disabled in analysis on AOD's");
       iPWG4PtQATPC        = 0;
+      if( iPWG4Cosmics)::Info("AnalysisTrainPWG4Jets.C::CheckModuleFlags", "PWG4 Comics disabled in analysis on AOD's");
+      iPWG4Cosmics        = 0;
+
       if( iPWG4PtSpectra)::Info("AnalysisTrainPWG4Jets.C::CheckModuleFlags", "PWG4 PtQAMC disabled in analysis on AOD's");
       iPWG4PtSpectra     = 0;
       if(iPWG4KMeans)::Info("AnalysisTrainPWG4Jets.C::CheckModuleFlags", "PWG4KMeans disabled on AOD's");
@@ -781,7 +796,7 @@ void CheckModuleFlags(const char *mode) {
        iPWG4JetSpectrum = iPWG4UE = iPWG4ThreeJets = iDIJETAN = 0;
       }
    }
-   iPWG4JetTasks = iPWG4JetServices||iPWG4JetSpectrum||iPWG4UE||iPWG4PtQAMC||iPWG4PtSpectra||iPWG4PtQATPC||iPWG4ThreeJets||iPWG4JetChem;
+   iPWG4JetTasks = iPWG4JetServices||iPWG4JetSpectrum||iPWG4UE||iPWG4PtQAMC||iPWG4PtSpectra||iPWG4PtQATPC||iPWG4Cosmics||iPWG4ThreeJets||iPWG4JetChem;
    iPWG4PartCorrLibs = iPWG4PartCorr||iPWG4Tagged||iPWG4CaloQA;
    iJETANLib = iPWG4JetTasks||iJETAN||iDIJETAN;
    if (iESDfilter) {iAODhandler=1;}