[u/mrichter/AliRoot.git] / PWGLF / spectraJET / AliAnalysisTaskJetSpectraAOD.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, 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.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//         AliAnalysisTaskJetSpectraAOD class
//-----------------------------------------------------------------

#include "TChain.h"
#include "TTree.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TH2F.h"
#include "THnSparse.h"
#include "TCanvas.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliVTrack.h"
#include "AliAODMCParticle.h"
#include "AliVParticle.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODInputHandler.h"
#include "AliAnalysisTaskJetSpectraAOD.h"
#include "AliAnalysisTaskESDfilter.h"
#include "AliAnalysisDataContainer.h"
#include "AliCentrality.h"
#include "TProof.h"
#include "AliVEvent.h"
#include "AliStack.h"
#include <TMCProcess.h>

#include "AliSpectraAODTrackCuts.h"
#include "AliSpectraAODEventCuts.h"

//jet include
#include "AliAODHandler.h"
#include "AliAODJetEventBackground.h"
#include "AliAODJet.h"
#include "AliAnalysisHelperJetTasks.h"

#include <iostream>

using namespace std;

ClassImp(AliAnalysisTaskJetSpectraAOD) 

//________________________________________________________________________
AliAnalysisTaskJetSpectraAOD::AliAnalysisTaskJetSpectraAOD(const char *name) : AliAnalysisTaskSE(name),
  fAOD(0),
  fIsMC(0),
  fEventCuts(0),
  fTrackCuts(0),
  fVZEROside(0),
  fOutput(0),
  fAODJets(0),
  fJetBranchName(""),
  fListJets(0),
  fBackgroundBranch(""),
  fOfflineTrgMask(AliVEvent::kMB),
  fFilterMask(0),
  fJetPtMin(0x0),
  fJetEtaMin(0x0),
  fJetEtaMax(0x0),
  fLeadPtMin(0x0),
  fnCentBins(20),
  fnQvecBins(20),
  fnptLeadBins(4),
  fIsQvecCalibMode(0),
  fQvecUpperLim(100),
  fIsQvecCut(0),
  fQvecMin(0),
  fQvecMax(100),
  fHistEvtSelection(0x0),
  fDebug(0),
  fMinNcontributors(0),
  fRejectPileup(0),
  fR(0.4),
  fZvertexDiff(1),
  fZvertex(10.)
{
  // Default constructor
  
  DefineInput(0, TChain::Class());
  DefineOutput(1, TList::Class());
  DefineOutput(2, AliSpectraAODEventCuts::Class());
  DefineOutput(3, AliSpectraAODTrackCuts::Class());
  
}

//________________________________________________________________________
AliAnalysisTaskJetSpectraAOD::~AliAnalysisTaskJetSpectraAOD()
{
   delete fListJets;
}
//________________________________________________________________________
//________________________________________________________________________
void AliAnalysisTaskJetSpectraAOD::UserCreateOutputObjects()
{
  // create output objects
  fOutput = new TList();
  fOutput->SetOwner();
  fOutput->SetName("chistpt");
  
  fListJets = new TList;
  
  if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
  if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
    
  // binning common to all the THn
  const Double_t ptBins[] = {-50.,-45.,-40.,-35.,-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,100.,150.,200.};
  const Int_t nptBins=31;
  
  //dimensions of THnSparse for jets
  const Int_t nvarjet=5;
  //                                        pt_raw    pt_corr           cent             Q vec                pt_lead
  Int_t    binsHistRealJet[nvarjet] = {    nptBins,   nptBins,       fnCentBins,      fnQvecBins,       fnptLeadBins};
  Double_t xminHistRealJet[nvarjet] = {         0.,        0.,             0.,                0.,                  0.};
  Double_t xmaxHistRealJet[nvarjet] = {       200.,      200.,           100.,     fQvecUpperLim,                 20.};    
  THnSparseF* NSparseHistJet = new THnSparseF("NSparseHistJet","NSparseHistJet",nvarjet,binsHistRealJet,xminHistRealJet,xmaxHistRealJet);
  NSparseHistJet->GetAxis(0)->SetTitle("#it{p}_{T,raw}");
  NSparseHistJet->GetAxis(0)->SetName("pT_raw");
  NSparseHistJet->SetBinEdges(0,ptBins);
  NSparseHistJet->GetAxis(1)->SetTitle("#it{p}_{T,corr}");
  NSparseHistJet->GetAxis(1)->SetName("pT_corr");
  NSparseHistJet->SetBinEdges(1,ptBins);
  NSparseHistJet->GetAxis(2)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistJet->GetAxis(2)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistJet->GetAxis(3)->SetTitle("Q vec");
  NSparseHistJet->GetAxis(3)->SetName("Q_vec");
  NSparseHistJet->GetAxis(4)->SetTitle("#it{p}_{T,lead}");
  NSparseHistJet->GetAxis(4)->SetName("pT_lead");
  fOutput->Add(NSparseHistJet);
  
  //dimensions of THnSparse for the normalization
  const Int_t nvarev=3;
  //                                             cent         Q vec         rho
  Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,      100 };
  Double_t xminHistRealEv[nvarev] = {           0.,               0.,        0.};
  Double_t xmaxHistRealEv[nvarev] = {         100.,      fQvecUpperLim,    200.};
  THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
  NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
  NSparseHistEv->GetAxis(1)->SetName("Q_vec");
  NSparseHistEv->GetAxis(2)->SetTitle("rho");
  NSparseHistEv->GetAxis(2)->SetName("rho");
  fOutput->Add(NSparseHistEv);
  
//   TH1F* fHistTest = new TH1F("fHistTest","fHistTest",nptBins-1,ptBins);
//   fOutput->Add(fHistTest);

  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 7.5);
  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"centrality (rejected)");
  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"multiplicity (rejected)");
  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"ESE (rejected)");
  fOutput->Add(fHistEvtSelection);
  
  PostData(1, fOutput  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  
}

//________________________________________________________________________
void AliAnalysisTaskJetSpectraAOD::UserExec(Option_t *)
{
  
  // check for jet branches
  if(!strlen(fJetBranchName.Data())){
    AliError("Jet branch name not set.");
    return;
  }
  
  // main event loop
  fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
  if (!fAOD) {
    AliWarning("ERROR: AliAODEvent not available \n");
    return;
  }
    
  if (strcmp(fAOD->ClassName(), "AliAODEvent"))
    {
      AliFatal("Not processing AODs");
    }
  
  TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
  if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
    fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
  }

  /* -- event selection -- */
  fHistEvtSelection->Fill(1); // number of events before event selection

  //jet service task event selection.
  Bool_t selected=kTRUE;
  selected = AliAnalysisHelperJetTasks::Selected();
  if(!selected){
    // no selection by the service task, we continue
    PostData(1,fOutput);
    PostData(2, fEventCuts);
    PostData(3, fTrackCuts);
  return;}
  
  // physics selection: this is now redundant, all should appear as accepted after service task selection
  AliInputEventHandler* inputHandler = (AliInputEventHandler*)
  ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
  //std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
  if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
    if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
    fHistEvtSelection->Fill(2);
    PostData(1,fOutput);
    PostData(2, fEventCuts);
    PostData(3, fTrackCuts);
    return;
  }
  
  AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
  Int_t nTracksPrim = primVtx->GetNContributors();
  
  if (fDebug) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
  if(nTracksPrim < fMinNcontributors){
    if (fDebug) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); 
    fHistEvtSelection->Fill(3);
    PostData(1,fOutput);
    PostData(2, fEventCuts);
    PostData(3, fTrackCuts);
    return;
  }
  
  TString primVtxName(primVtx->GetName());

  if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
    if (fDebug) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
    fHistEvtSelection->Fill(4);
    PostData(1,fOutput);
    PostData(2, fEventCuts);
    PostData(3, fTrackCuts);
    return;
  }

  if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
    if (fDebug) Printf("%s:%d SPD pileup: event REJECTED...",(char*)__FILE__,__LINE__);
    fHistEvtSelection->Fill(5);
    PostData(1,fOutput);
    PostData(2, fEventCuts);
    PostData(3, fTrackCuts);
    return;
  }
    
  //cut on difference of the z-vertex
  if (fZvertexDiff || (fZvertex>0)) {
    
    Double_t vzPRI = +999;
    Double_t vzSPD = -999;
    
    const AliVVertex *pv = fAOD->GetPrimaryVertex(); // AOD
    if (pv && pv->GetNContributors()>0) {
      vzPRI = pv->GetZ();
    }
    
    const AliVVertex *sv = fAOD->GetPrimaryVertexSPD();
    if (sv && sv->GetNContributors()>0) {
      vzSPD = sv->GetZ();
    }
    
    Double_t  dvertex = TMath::Abs(vzPRI-vzSPD);

    if (fZvertexDiff && (dvertex>0.1)) {
      if (fDebug) Printf("%s:%d ZvertexDiff Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
      return;}
    if ((fZvertex>0) && (TMath::Abs(vzPRI)>fZvertex)) {
      if (fDebug) Printf("%s:%d ZvertexDiff Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
      return;}
  }
  
  //ESE event cuts
  if(!fEventCuts->IsSelected(fAOD,fTrackCuts)){
    if (fDebug) Printf("%s:%d ESE Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
    fHistEvtSelection->Fill(6);
    PostData(1,fOutput);
    PostData(2, fEventCuts);
    PostData(3, fTrackCuts);
    return;
  }
  
  if (fDebug) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__); 
  fHistEvtSelection->Fill(0.);
  // accepted events  
  // -- end event selection --
  
  
  Double_t Qvec=0.;//in case of MC we save space in the memory
  if(!fIsMC){
    if(fIsQvecCalibMode){
      if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
      else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
    }
    else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
  }
  
  if(fIsQvecCut && (Qvec<fQvecMin || Qvec>fQvecMax) ) return;
  
  Double_t Cent=fEventCuts->GetCent();

  // get background
  AliAODJetEventBackground* externalBackground = 0;
  if(fAODJets && !externalBackground && fBackgroundBranch.Length()){
    externalBackground =  (AliAODJetEventBackground*)(fAODJets->FindListObject(fBackgroundBranch.Data()));
    if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
  }
  
  Float_t rho = 0;
  if(externalBackground)rho = externalBackground->GetBackground(0);  //default schema
    if(rho==0) rho=-9999.; //rho value = 0 are non-physical -> removed from the distribution
  
  // fetch jets 
  TClonesArray *aodJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fJetBranchName.Data()));
  if(!aodJets){
    AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName.Data()));
    if(fAOD){
      Printf("Input AOD >>>>");
      fAOD->Print();    
    }
    return;
  }
  
  Bool_t kDeltaPt = kFALSE;
  if(fJetBranchName.Contains("RandomCone")){
//     cout<<"RANDOM CONES BRANCH!"<<endl;
    kDeltaPt = kTRUE;
  }
  
  
  fListJets->Clear();
  for (Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
    AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
    if (jet) fListJets->Add(jet);
  }
  

  for(Int_t i=0; i<fListJets->GetEntries(); ++i){
    AliAODJet* jet = (AliAODJet*)(fListJets->At(i));
    
    if((jet->Eta()<fJetEtaMin)||(jet->Eta()>fJetEtaMax)) continue;
    
    Double_t ptJet   = jet->Pt();
    
    Double_t areaJet = jet->EffectiveAreaCharged();

    Double_t ptcorr = -999.;
    if(kDeltaPt) ptcorr = ptJet - (rho*TMath::Pi()*fR*fR);
    else ptcorr = ptJet-(rho*areaJet);
           
    Double_t varJet[5];
    varJet[0]=jet->Pt();
    varJet[1]=(Double_t)ptcorr;
    varJet[2]=(Double_t)Cent;
    varJet[3]=(Double_t)Qvec;
    
    if(!kDeltaPt){
      AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
      if(fLeadPtMin>0 && leadTrack->Pt()<fLeadPtMin)continue;
    
      varJet[4]=(Double_t)leadTrack->Pt();
    }
    else varJet[4] = -1;
    
    ((THnSparseF*)fOutput->FindObject("NSparseHistJet"))->Fill(varJet);//jet loop
  }
  
  
//   //track loop
//   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {  //FIXME loop on aod track... should be on jet tracks???
//     AliAODTrack* track = fAOD->GetTrack(iTracks);
// //     if(!(track->TestFilterBit(1024)))continue;
//     if (!fTrackCuts->IsSelected(track,kTRUE)) continue;
//       
//     TH1F* h=(TH1F*)fOutput->FindObject("fHistTest");h->Fill(track->Pt());
//     
//   } // end loop on tracks
  
  Double_t varEv[3];
  varEv[0]=Cent;
  varEv[1]=Qvec;
  varEv[2]=rho;
  ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
  
  PostData(1,fOutput);
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  //Printf("............. end of Exec");
  
}

//_________________________________________________________________
AliVParticle *AliAnalysisTaskJetSpectraAOD::LeadingTrackFromJetRefs(AliAODJet* jet){
  if(!jet)return 0;
  TRefArray *refs = jet->GetRefTracks();
  if(!refs) return 0;
  AliVParticle *leading = 0;
  Float_t fMaxPt = 0;
  for(int i = 0;i<refs->GetEntriesFast();i++){
    AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
    if(!tmp)continue;
    if(tmp->Pt()>fMaxPt){
      leading = tmp;
      fMaxPt = tmp->Pt();
    }
  }
  return leading;
}

//_________________________________________________________________
void   AliAnalysisTaskJetSpectraAOD::Terminate(Option_t *)
{
  // Terminate
}

//jet
void AliAnalysisTaskJetSpectraAOD::SetBranchNames(const TString &branch)
{
   fJetBranchName = branch;
}

void AliAnalysisTaskJetSpectraAOD::SetRecBackgroundBranch(const TString &bckbranch)
{
   fBackgroundBranch = bckbranch;
}
Commit	Line	Data
582efd60	1	/**************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-----------------------------------------------------------------
	17	// AliAnalysisTaskJetSpectraAOD class
	18	//-----------------------------------------------------------------
	19
	20	#include "TChain.h"
	21	#include "TTree.h"
	22	#include "TLegend.h"
	23	#include "TH1F.h"
	24	#include "TH2F.h"
	25	#include "THnSparse.h"
	26	#include "TCanvas.h"
	27	#include "AliAnalysisTask.h"
	28	#include "AliAnalysisManager.h"
	29	#include "AliVTrack.h"
	30	#include "AliAODMCParticle.h"
	31	#include "AliVParticle.h"
	32	#include "AliAODEvent.h"
	33	#include "AliAODTrack.h"
	34	#include "AliAODInputHandler.h"
	35	#include "AliAnalysisTaskJetSpectraAOD.h"
	36	#include "AliAnalysisTaskESDfilter.h"
	37	#include "AliAnalysisDataContainer.h"
	38	#include "AliCentrality.h"
	39	#include "TProof.h"
	40	#include "AliVEvent.h"
	41	#include "AliStack.h"
	42	#include <TMCProcess.h>
	43
	44	#include "AliSpectraAODTrackCuts.h"
	45	#include "AliSpectraAODEventCuts.h"
	46
	47	//jet include
	48	#include "AliAODHandler.h"
	49	#include "AliAODJetEventBackground.h"
	50	#include "AliAODJet.h"
296131cf	51	#include "AliAnalysisHelperJetTasks.h"
582efd60	52
	53	#include <iostream>
	54
	55	using namespace std;
	56
	57	ClassImp(AliAnalysisTaskJetSpectraAOD)
	58
	59	//________________________________________________________________________
	60	AliAnalysisTaskJetSpectraAOD::AliAnalysisTaskJetSpectraAOD(const char *name) : AliAnalysisTaskSE(name),
	61	fAOD(0),
	62	fIsMC(0),
	63	fEventCuts(0),
	64	fTrackCuts(0),
	65	fVZEROside(0),
	66	fOutput(0),
	67	fAODJets(0),
	68	fJetBranchName(""),
	69	fListJets(0),
	70	fBackgroundBranch(""),
296131cf	71	fOfflineTrgMask(AliVEvent::kMB),
582efd60	72	fFilterMask(0),
296131cf	73	fJetPtMin(0x0),
582efd60	74	fJetEtaMin(0x0),
582efd60	75	fJetEtaMax(0x0),
296131cf	76	fLeadPtMin(0x0),
582efd60	77	fnCentBins(20),
582efd60	78	fnQvecBins(20),
296131cf	79	fnptLeadBins(4),
582efd60	80	fIsQvecCalibMode(0),
	81	fQvecUpperLim(100),
	82	fIsQvecCut(0),
	83	fQvecMin(0),
296131cf	84	fQvecMax(100),
	85	fHistEvtSelection(0x0),
	86	fDebug(0),
	87	fMinNcontributors(0),
9a924571	88	fRejectPileup(0),
1b8d95a7	89	fR(0.4),
	90	fZvertexDiff(1),
	91	fZvertex(10.)
582efd60	92	{
	93	// Default constructor
	94
	95	DefineInput(0, TChain::Class());
	96	DefineOutput(1, TList::Class());
	97	DefineOutput(2, AliSpectraAODEventCuts::Class());
	98	DefineOutput(3, AliSpectraAODTrackCuts::Class());
	99
	100	}
	101
	102	//________________________________________________________________________
	103	AliAnalysisTaskJetSpectraAOD::~AliAnalysisTaskJetSpectraAOD()
	104	{
	105	delete fListJets;
	106	}
	107	//________________________________________________________________________
	108	//________________________________________________________________________
	109	void AliAnalysisTaskJetSpectraAOD::UserCreateOutputObjects()
	110	{
9a924571	111	// create output objects
582efd60	112	fOutput = new TList();
	113	fOutput->SetOwner();
	114	fOutput->SetName("chistpt");
	115
	116	fListJets = new TList;
	117
	118	if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
	119	if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
	120
	121	// binning common to all the THn
9a924571	122	const Double_t ptBins[] = {-50.,-45.,-40.,-35.,-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,100.,150.,200.};
9a924571	123	const Int_t nptBins=31;
582efd60	124
582efd60	125	//dimensions of THnSparse for jets
9a924571	126	const Int_t nvarjet=5;
	127	// pt_raw pt_corr cent Q vec pt_lead
	128	Int_t binsHistRealJet[nvarjet] = { nptBins, nptBins, fnCentBins, fnQvecBins, fnptLeadBins};
	129	Double_t xminHistRealJet[nvarjet] = { 0., 0., 0., 0., 0.};
	130	Double_t xmaxHistRealJet[nvarjet] = { 200., 200., 100., fQvecUpperLim, 20.};
582efd60	131	THnSparseF* NSparseHistJet = new THnSparseF("NSparseHistJet","NSparseHistJet",nvarjet,binsHistRealJet,xminHistRealJet,xmaxHistRealJet);
	132	NSparseHistJet->GetAxis(0)->SetTitle("#it{p}_{T,raw}");
	133	NSparseHistJet->GetAxis(0)->SetName("pT_raw");
	134	NSparseHistJet->SetBinEdges(0,ptBins);
	135	NSparseHistJet->GetAxis(1)->SetTitle("#it{p}_{T,corr}");
	136	NSparseHistJet->GetAxis(1)->SetName("pT_corr");
9a924571	137	NSparseHistJet->SetBinEdges(1,ptBins);
582efd60	138	NSparseHistJet->GetAxis(2)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	139	NSparseHistJet->GetAxis(2)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	140	NSparseHistJet->GetAxis(3)->SetTitle("Q vec");
	141	NSparseHistJet->GetAxis(3)->SetName("Q_vec");
9a924571	142	NSparseHistJet->GetAxis(4)->SetTitle("#it{p}_{T,lead}");
9a924571	143	NSparseHistJet->GetAxis(4)->SetName("pT_lead");
582efd60	144	fOutput->Add(NSparseHistJet);
	145
	146	//dimensions of THnSparse for the normalization
	147	const Int_t nvarev=3;
	148	// cent Q vec rho
2ee86b44	149	Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, 100 };
582efd60	150	Double_t xminHistRealEv[nvarev] = { 0., 0., 0.};
	151	Double_t xmaxHistRealEv[nvarev] = { 100., fQvecUpperLim, 200.};
	152	THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
	153	NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	154	NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	155	NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
	156	NSparseHistEv->GetAxis(1)->SetName("Q_vec");
	157	NSparseHistEv->GetAxis(2)->SetTitle("rho");
	158	NSparseHistEv->GetAxis(2)->SetName("rho");
	159	fOutput->Add(NSparseHistEv);
	160
	161	// TH1F* fHistTest = new TH1F("fHistTest","fHistTest",nptBins-1,ptBins);
	162	// fOutput->Add(fHistTest);
296131cf	163
	164	fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 7.5);
	165	fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
	166	fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
	167	fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
	168	fHistEvtSelection->GetXaxis()->SetBinLabel(4,"centrality (rejected)");
	169	fHistEvtSelection->GetXaxis()->SetBinLabel(5,"multiplicity (rejected)");
	170	fHistEvtSelection->GetXaxis()->SetBinLabel(6,"ESE (rejected)");
	171	fOutput->Add(fHistEvtSelection);
582efd60	172
	173	PostData(1, fOutput );
	174	PostData(2, fEventCuts);
	175	PostData(3, fTrackCuts);
	176
	177	}
	178
	179	//________________________________________________________________________
	180	void AliAnalysisTaskJetSpectraAOD::UserExec(Option_t *)
	181	{
	182
	183	// check for jet branches
	184	if(!strlen(fJetBranchName.Data())){
	185	AliError("Jet branch name not set.");
	186	return;
	187	}
	188
	189	// main event loop
	190	fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
	191	if (!fAOD) {
	192	AliWarning("ERROR: AliAODEvent not available \n");
	193	return;
	194	}
	195
	196	if (strcmp(fAOD->ClassName(), "AliAODEvent"))
	197	{
	198	AliFatal("Not processing AODs");
	199	}
	200
	201	TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
	202	if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
	203	fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
	204	}
	205
296131cf	206	/* -- event selection -- */
296131cf	207	fHistEvtSelection->Fill(1); // number of events before event selection
582efd60	208
296131cf	209	//jet service task event selection.
	210	Bool_t selected=kTRUE;
	211	selected = AliAnalysisHelperJetTasks::Selected();
	212	if(!selected){
	213	// no selection by the service task, we continue
	214	PostData(1,fOutput);
	215	PostData(2, fEventCuts);
	216	PostData(3, fTrackCuts);
	217	return;}
	218
	219	// physics selection: this is now redundant, all should appear as accepted after service task selection
	220	AliInputEventHandler* inputHandler = (AliInputEventHandler*)
	221	((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
	222	//std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
	223	if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
	224	if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
	225	fHistEvtSelection->Fill(2);
	226	PostData(1,fOutput);
	227	PostData(2, fEventCuts);
	228	PostData(3, fTrackCuts);
	229	return;
	230	}
	231
	232	AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
	233	Int_t nTracksPrim = primVtx->GetNContributors();
	234
	235	if (fDebug) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
	236	if(nTracksPrim < fMinNcontributors){
	237	if (fDebug) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
	238	fHistEvtSelection->Fill(3);
	239	PostData(1,fOutput);
	240	PostData(2, fEventCuts);
	241	PostData(3, fTrackCuts);
	242	return;
	243	}
	244
	245	TString primVtxName(primVtx->GetName());
	246
	247	if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
	248	if (fDebug) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
	249	fHistEvtSelection->Fill(4);
	250	PostData(1,fOutput);
	251	PostData(2, fEventCuts);
	252	PostData(3, fTrackCuts);
	253	return;
	254	}
	255
	256	if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
	257	if (fDebug) Printf("%s:%d SPD pileup: event REJECTED...",(char*)__FILE__,__LINE__);
	258	fHistEvtSelection->Fill(5);
	259	PostData(1,fOutput);
	260	PostData(2, fEventCuts);
	261	PostData(3, fTrackCuts);
	262	return;
	263	}
1b8d95a7	264
	265	//cut on difference of the z-vertex
	266	if (fZvertexDiff \|\| (fZvertex>0)) {
	267
	268	Double_t vzPRI = +999;
	269	Double_t vzSPD = -999;
	270
	271	const AliVVertex *pv = fAOD->GetPrimaryVertex(); // AOD
	272	if (pv && pv->GetNContributors()>0) {
	273	vzPRI = pv->GetZ();
	274	}
	275
	276	const AliVVertex *sv = fAOD->GetPrimaryVertexSPD();
	277	if (sv && sv->GetNContributors()>0) {
	278	vzSPD = sv->GetZ();
	279	}
	280
	281	Double_t dvertex = TMath::Abs(vzPRI-vzSPD);
296131cf	282
1b8d95a7	283	if (fZvertexDiff && (dvertex>0.1)) {
	284	if (fDebug) Printf("%s:%d ZvertexDiff Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
	285	return;}
	286	if ((fZvertex>0) && (TMath::Abs(vzPRI)>fZvertex)) {
	287	if (fDebug) Printf("%s:%d ZvertexDiff Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
	288	return;}
	289	}
	290
296131cf	291	//ESE event cuts
	292	if(!fEventCuts->IsSelected(fAOD,fTrackCuts)){
	293	if (fDebug) Printf("%s:%d ESE Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
	294	fHistEvtSelection->Fill(6);
	295	PostData(1,fOutput);
	296	PostData(2, fEventCuts);
	297	PostData(3, fTrackCuts);
	298	return;
	299	}
	300
	301	if (fDebug) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
	302	fHistEvtSelection->Fill(0.);
	303	// accepted events
	304	// -- end event selection --
	305
582efd60	306
	307	Double_t Qvec=0.;//in case of MC we save space in the memory
	308	if(!fIsMC){
	309	if(fIsQvecCalibMode){
	310	if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
	311	else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
	312	}
	313	else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
	314	}
	315
	316	if(fIsQvecCut && (Qvec<fQvecMin \|\| Qvec>fQvecMax) ) return;
	317
	318	Double_t Cent=fEventCuts->GetCent();
296131cf	319
582efd60	320	// get background
	321	AliAODJetEventBackground* externalBackground = 0;
	322	if(fAODJets && !externalBackground && fBackgroundBranch.Length()){
	323	externalBackground = (AliAODJetEventBackground*)(fAODJets->FindListObject(fBackgroundBranch.Data()));
	324	if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
	325	}
	326
	327	Float_t rho = 0;
	328	if(externalBackground)rho = externalBackground->GetBackground(0); //default schema
72ffccdd	329	if(rho==0) rho=-9999.; //rho value = 0 are non-physical -> removed from the distribution
582efd60	330
2ee86b44	331	// fetch jets
582efd60	332	TClonesArray aodJets = dynamic_cast<TClonesArray>(fAODJets->FindListObject(fJetBranchName.Data()));
	333	if(!aodJets){
	334	AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName.Data()));
	335	if(fAOD){
	336	Printf("Input AOD >>>>");
	337	fAOD->Print();
	338	}
	339	return;
	340	}
	341
9a924571	342	Bool_t kDeltaPt = kFALSE;
	343	if(fJetBranchName.Contains("RandomCone")){
	344	// cout<<"RANDOM CONES BRANCH!"<<endl;
	345	kDeltaPt = kTRUE;
	346	}
	347
	348
582efd60	349	fListJets->Clear();
	350	for (Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
	351	AliAODJet jet = dynamic_cast<AliAODJet>((*aodJets)[iJet]);
	352	if (jet) fListJets->Add(jet);
	353	}
	354
	355
	356
	357	for(Int_t i=0; i<fListJets->GetEntries(); ++i){
	358	AliAODJet* jet = (AliAODJet*)(fListJets->At(i));
	359
	360	if((jet->Eta()<fJetEtaMin)\|\|(jet->Eta()>fJetEtaMax)) continue;
	361
	362	Double_t ptJet = jet->Pt();
	363
	364	Double_t areaJet = jet->EffectiveAreaCharged();
	365
9a924571	366	Double_t ptcorr = -999.;
	367	if(kDeltaPt) ptcorr = ptJet - (rhoTMath::Pi()fR*fR);
	368	else ptcorr = ptJet-(rho*areaJet);
582efd60	369
9a924571	370	Double_t varJet[5];
582efd60	371	varJet[0]=jet->Pt();
	372	varJet[1]=(Double_t)ptcorr;
	373	varJet[2]=(Double_t)Cent;
	374	varJet[3]=(Double_t)Qvec;
582efd60	375
9a924571	376	if(!kDeltaPt){
	377	AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
	378	if(fLeadPtMin>0 && leadTrack->Pt()<fLeadPtMin)continue;
296131cf	379
9a924571	380	varJet[4]=(Double_t)leadTrack->Pt();
	381	}
	382	else varJet[4] = -1;
296131cf	383
582efd60	384	((THnSparseF*)fOutput->FindObject("NSparseHistJet"))->Fill(varJet);//jet loop
	385	}
	386
	387
	388	// //track loop
	389	// for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) { //FIXME loop on aod track... should be on jet tracks???
	390	// AliAODTrack* track = fAOD->GetTrack(iTracks);
	391	// // if(!(track->TestFilterBit(1024)))continue;
	392	// if (!fTrackCuts->IsSelected(track,kTRUE)) continue;
	393	//
	394	// TH1F* h=(TH1F*)fOutput->FindObject("fHistTest");h->Fill(track->Pt());
	395	//
	396	// } // end loop on tracks
	397
	398	Double_t varEv[3];
	399	varEv[0]=Cent;
	400	varEv[1]=Qvec;
	401	varEv[2]=rho;
	402	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
	403
	404	PostData(1,fOutput);
	405	PostData(2, fEventCuts);
	406	PostData(3, fTrackCuts);
	407	//Printf("............. end of Exec");
	408
	409	}
296131cf	410
	411	//_________________________________________________________________
	412	AliVParticle AliAnalysisTaskJetSpectraAOD::LeadingTrackFromJetRefs(AliAODJet jet){
	413	if(!jet)return 0;
	414	TRefArray *refs = jet->GetRefTracks();
	415	if(!refs) return 0;
	416	AliVParticle *leading = 0;
	417	Float_t fMaxPt = 0;
	418	for(int i = 0;i<refs->GetEntriesFast();i++){
	419	AliVParticle tmp = dynamic_cast<AliVParticle>(refs->At(i));
	420	if(!tmp)continue;
	421	if(tmp->Pt()>fMaxPt){
	422	leading = tmp;
	423	fMaxPt = tmp->Pt();
	424	}
	425	}
	426	return leading;
	427	}
	428
582efd60	429	//_________________________________________________________________
	430	void AliAnalysisTaskJetSpectraAOD::Terminate(Option_t *)
	431	{
	432	// Terminate
582efd60	433	}
	434
	435	//jet
	436	void AliAnalysisTaskJetSpectraAOD::SetBranchNames(const TString &branch)
	437	{
	438	fJetBranchName = branch;
	439	}
	440
	441	void AliAnalysisTaskJetSpectraAOD::SetRecBackgroundBranch(const TString &bckbranch)
	442	{
	443	fBackgroundBranch = bckbranch;
296131cf	444	}