[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / 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)
{
  // 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;
  }

  //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),
9a924571	89	fR(0.4)
582efd60	90	{
	91	// Default constructor
	92
	93	DefineInput(0, TChain::Class());
	94	DefineOutput(1, TList::Class());
	95	DefineOutput(2, AliSpectraAODEventCuts::Class());
	96	DefineOutput(3, AliSpectraAODTrackCuts::Class());
	97
	98	}
	99
	100	//________________________________________________________________________
	101	AliAnalysisTaskJetSpectraAOD::~AliAnalysisTaskJetSpectraAOD()
	102	{
	103	delete fListJets;
	104	}
	105	//________________________________________________________________________
	106	//________________________________________________________________________
	107	void AliAnalysisTaskJetSpectraAOD::UserCreateOutputObjects()
	108	{
9a924571	109	// create output objects
582efd60	110	fOutput = new TList();
	111	fOutput->SetOwner();
	112	fOutput->SetName("chistpt");
	113
	114	fListJets = new TList;
	115
	116	if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
	117	if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
	118
	119	// binning common to all the THn
9a924571	120	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	121	const Int_t nptBins=31;
582efd60	122
582efd60	123	//dimensions of THnSparse for jets
9a924571	124	const Int_t nvarjet=5;
	125	// pt_raw pt_corr cent Q vec pt_lead
	126	Int_t binsHistRealJet[nvarjet] = { nptBins, nptBins, fnCentBins, fnQvecBins, fnptLeadBins};
	127	Double_t xminHistRealJet[nvarjet] = { 0., 0., 0., 0., 0.};
	128	Double_t xmaxHistRealJet[nvarjet] = { 200., 200., 100., fQvecUpperLim, 20.};
582efd60	129	THnSparseF* NSparseHistJet = new THnSparseF("NSparseHistJet","NSparseHistJet",nvarjet,binsHistRealJet,xminHistRealJet,xmaxHistRealJet);
	130	NSparseHistJet->GetAxis(0)->SetTitle("#it{p}_{T,raw}");
	131	NSparseHistJet->GetAxis(0)->SetName("pT_raw");
	132	NSparseHistJet->SetBinEdges(0,ptBins);
	133	NSparseHistJet->GetAxis(1)->SetTitle("#it{p}_{T,corr}");
	134	NSparseHistJet->GetAxis(1)->SetName("pT_corr");
9a924571	135	NSparseHistJet->SetBinEdges(1,ptBins);
582efd60	136	NSparseHistJet->GetAxis(2)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	137	NSparseHistJet->GetAxis(2)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	138	NSparseHistJet->GetAxis(3)->SetTitle("Q vec");
	139	NSparseHistJet->GetAxis(3)->SetName("Q_vec");
9a924571	140	NSparseHistJet->GetAxis(4)->SetTitle("#it{p}_{T,lead}");
9a924571	141	NSparseHistJet->GetAxis(4)->SetName("pT_lead");
582efd60	142	fOutput->Add(NSparseHistJet);
	143
	144	//dimensions of THnSparse for the normalization
	145	const Int_t nvarev=3;
	146	// cent Q vec rho
2ee86b44	147	Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, 100 };
582efd60	148	Double_t xminHistRealEv[nvarev] = { 0., 0., 0.};
	149	Double_t xmaxHistRealEv[nvarev] = { 100., fQvecUpperLim, 200.};
	150	THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
	151	NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	152	NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	153	NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
	154	NSparseHistEv->GetAxis(1)->SetName("Q_vec");
	155	NSparseHistEv->GetAxis(2)->SetTitle("rho");
	156	NSparseHistEv->GetAxis(2)->SetName("rho");
	157	fOutput->Add(NSparseHistEv);
	158
	159	// TH1F* fHistTest = new TH1F("fHistTest","fHistTest",nptBins-1,ptBins);
	160	// fOutput->Add(fHistTest);
296131cf	161
	162	fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 7.5);
	163	fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
	164	fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
	165	fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
	166	fHistEvtSelection->GetXaxis()->SetBinLabel(4,"centrality (rejected)");
	167	fHistEvtSelection->GetXaxis()->SetBinLabel(5,"multiplicity (rejected)");
	168	fHistEvtSelection->GetXaxis()->SetBinLabel(6,"ESE (rejected)");
	169	fOutput->Add(fHistEvtSelection);
582efd60	170
	171	PostData(1, fOutput );
	172	PostData(2, fEventCuts);
	173	PostData(3, fTrackCuts);
	174
	175	}
	176
	177	//________________________________________________________________________
	178	void AliAnalysisTaskJetSpectraAOD::UserExec(Option_t *)
	179	{
	180
	181	// check for jet branches
	182	if(!strlen(fJetBranchName.Data())){
	183	AliError("Jet branch name not set.");
	184	return;
	185	}
	186
	187	// main event loop
	188	fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
	189	if (!fAOD) {
	190	AliWarning("ERROR: AliAODEvent not available \n");
	191	return;
	192	}
	193
	194	if (strcmp(fAOD->ClassName(), "AliAODEvent"))
	195	{
	196	AliFatal("Not processing AODs");
	197	}
	198
	199	TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
	200	if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
	201	fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
	202	}
	203
296131cf	204	/* -- event selection -- */
296131cf	205	fHistEvtSelection->Fill(1); // number of events before event selection
582efd60	206
296131cf	207	//jet service task event selection.
	208	Bool_t selected=kTRUE;
	209	selected = AliAnalysisHelperJetTasks::Selected();
	210	if(!selected){
	211	// no selection by the service task, we continue
	212	PostData(1,fOutput);
	213	PostData(2, fEventCuts);
	214	PostData(3, fTrackCuts);
	215	return;}
	216
	217	// physics selection: this is now redundant, all should appear as accepted after service task selection
	218	AliInputEventHandler* inputHandler = (AliInputEventHandler*)
	219	((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
	220	//std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
	221	if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
	222	if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
	223	fHistEvtSelection->Fill(2);
	224	PostData(1,fOutput);
	225	PostData(2, fEventCuts);
	226	PostData(3, fTrackCuts);
	227	return;
	228	}
	229
	230	AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
	231	Int_t nTracksPrim = primVtx->GetNContributors();
	232
	233	if (fDebug) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
	234	if(nTracksPrim < fMinNcontributors){
	235	if (fDebug) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
	236	fHistEvtSelection->Fill(3);
	237	PostData(1,fOutput);
	238	PostData(2, fEventCuts);
	239	PostData(3, fTrackCuts);
	240	return;
	241	}
	242
	243	TString primVtxName(primVtx->GetName());
	244
	245	if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
	246	if (fDebug) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
	247	fHistEvtSelection->Fill(4);
	248	PostData(1,fOutput);
	249	PostData(2, fEventCuts);
	250	PostData(3, fTrackCuts);
	251	return;
	252	}
	253
	254	if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
	255	if (fDebug) Printf("%s:%d SPD pileup: event REJECTED...",(char*)__FILE__,__LINE__);
	256	fHistEvtSelection->Fill(5);
	257	PostData(1,fOutput);
	258	PostData(2, fEventCuts);
	259	PostData(3, fTrackCuts);
	260	return;
	261	}
	262
	263	//ESE event cuts
	264	if(!fEventCuts->IsSelected(fAOD,fTrackCuts)){
	265	if (fDebug) Printf("%s:%d ESE Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
	266	fHistEvtSelection->Fill(6);
	267	PostData(1,fOutput);
	268	PostData(2, fEventCuts);
	269	PostData(3, fTrackCuts);
	270	return;
271	}
272
273	if (fDebug) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
274	fHistEvtSelection->Fill(0.);
275	// accepted events
276	// -- end event selection --
277
582efd60	278
	279	Double_t Qvec=0.;//in case of MC we save space in the memory
	280	if(!fIsMC){
	281	if(fIsQvecCalibMode){
	282	if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
	283	else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
	284	}
	285	else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
	286	}
	287
	288	if(fIsQvecCut && (Qvec<fQvecMin \|\| Qvec>fQvecMax) ) return;
	289
	290	Double_t Cent=fEventCuts->GetCent();
296131cf	291
582efd60	292	// get background
	293	AliAODJetEventBackground* externalBackground = 0;
	294	if(fAODJets && !externalBackground && fBackgroundBranch.Length()){
	295	externalBackground = (AliAODJetEventBackground*)(fAODJets->FindListObject(fBackgroundBranch.Data()));
	296	if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
	297	}
	298
	299	Float_t rho = 0;
	300	if(externalBackground)rho = externalBackground->GetBackground(0); //default schema
72ffccdd	301	if(rho==0) rho=-9999.; //rho value = 0 are non-physical -> removed from the distribution
582efd60	302
2ee86b44	303	// fetch jets
582efd60	304	TClonesArray aodJets = dynamic_cast<TClonesArray>(fAODJets->FindListObject(fJetBranchName.Data()));
	305	if(!aodJets){
	306	AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName.Data()));
	307	if(fAOD){
	308	Printf("Input AOD >>>>");
	309	fAOD->Print();
	310	}
	311	return;
	312	}
	313
9a924571	314	Bool_t kDeltaPt = kFALSE;
	315	if(fJetBranchName.Contains("RandomCone")){
	316	// cout<<"RANDOM CONES BRANCH!"<<endl;
	317	kDeltaPt = kTRUE;
	318	}
	319
	320
582efd60	321	fListJets->Clear();
	322	for (Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
	323	AliAODJet jet = dynamic_cast<AliAODJet>((*aodJets)[iJet]);
	324	if (jet) fListJets->Add(jet);
	325	}
	326
	327
	328
	329	for(Int_t i=0; i<fListJets->GetEntries(); ++i){
	330	AliAODJet* jet = (AliAODJet*)(fListJets->At(i));
	331
	332	if((jet->Eta()<fJetEtaMin)\|\|(jet->Eta()>fJetEtaMax)) continue;
	333
	334	Double_t ptJet = jet->Pt();
	335
	336	Double_t areaJet = jet->EffectiveAreaCharged();
	337
9a924571	338	Double_t ptcorr = -999.;
	339	if(kDeltaPt) ptcorr = ptJet - (rhoTMath::Pi()fR*fR);
	340	else ptcorr = ptJet-(rho*areaJet);
582efd60	341
9a924571	342	Double_t varJet[5];
582efd60	343	varJet[0]=jet->Pt();
	344	varJet[1]=(Double_t)ptcorr;
	345	varJet[2]=(Double_t)Cent;
	346	varJet[3]=(Double_t)Qvec;
582efd60	347
9a924571	348	if(!kDeltaPt){
	349	AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
	350	if(fLeadPtMin>0 && leadTrack->Pt()<fLeadPtMin)continue;
296131cf	351
9a924571	352	varJet[4]=(Double_t)leadTrack->Pt();
	353	}
	354	else varJet[4] = -1;
296131cf	355
582efd60	356	((THnSparseF*)fOutput->FindObject("NSparseHistJet"))->Fill(varJet);//jet loop
	357	}
	358
	359
	360	// //track loop
	361	// for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) { //FIXME loop on aod track... should be on jet tracks???
	362	// AliAODTrack* track = fAOD->GetTrack(iTracks);
	363	// // if(!(track->TestFilterBit(1024)))continue;
	364	// if (!fTrackCuts->IsSelected(track,kTRUE)) continue;
	365	//
	366	// TH1F* h=(TH1F*)fOutput->FindObject("fHistTest");h->Fill(track->Pt());
	367	//
	368	// } // end loop on tracks
	369
	370	Double_t varEv[3];
	371	varEv[0]=Cent;
	372	varEv[1]=Qvec;
	373	varEv[2]=rho;
	374	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
	375
	376	PostData(1,fOutput);
	377	PostData(2, fEventCuts);
	378	PostData(3, fTrackCuts);
	379	//Printf("............. end of Exec");
	380
	381	}
296131cf	382
	383	//_________________________________________________________________
	384	AliVParticle AliAnalysisTaskJetSpectraAOD::LeadingTrackFromJetRefs(AliAODJet jet){
	385	if(!jet)return 0;
	386	TRefArray *refs = jet->GetRefTracks();
	387	if(!refs) return 0;
	388	AliVParticle *leading = 0;
	389	Float_t fMaxPt = 0;
	390	for(int i = 0;i<refs->GetEntriesFast();i++){
	391	AliVParticle tmp = dynamic_cast<AliVParticle>(refs->At(i));
	392	if(!tmp)continue;
	393	if(tmp->Pt()>fMaxPt){
	394	leading = tmp;
	395	fMaxPt = tmp->Pt();
	396	}
	397	}
	398	return leading;
	399	}
	400
582efd60	401	//_________________________________________________________________
	402	void AliAnalysisTaskJetSpectraAOD::Terminate(Option_t *)
	403	{
	404	// Terminate
582efd60	405	}
	406
	407	//jet
	408	void AliAnalysisTaskJetSpectraAOD::SetBranchNames(const TString &branch)
	409	{
	410	fJetBranchName = branch;
	411	}
	412
	413	void AliAnalysisTaskJetSpectraAOD::SetRecBackgroundBranch(const TString &bckbranch)
	414	{
	415	fBackgroundBranch = bckbranch;
296131cf	416	}