[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)
{
  // 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()
{
  Printf("\n\n\n\n\n\n In CreateOutput Object:");
  
  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[] = {0.15,5.,10.,15.,20.,25.,30.,35.,40.,50.,75.,100.,150.,200.};
  const Int_t nptBins=13;
  
  //dimensions of THnSparse for jets
  const Int_t nvarjet=6;
  //                                        pt_raw    pt_corr           cent             Q vec            rho          pt_lead
  Int_t    binsHistRealJet[nvarjet] = {    nptBins,   nptBins,       fnCentBins,      fnQvecBins,          40.,    fnptLeadBins};
  Double_t xminHistRealJet[nvarjet] = {         0.,        0.,             0.,                0.,           0.,              0.};
  Double_t xmaxHistRealJet[nvarjet] = {       200.,      200.,           100.,     fQvecUpperLim,         200.,             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("rho");
  NSparseHistJet->GetAxis(4)->SetName("rho");
  NSparseHistJet->GetAxis(5)->SetTitle("#it{p}_{T,lead}");
  NSparseHistJet->GetAxis(5)->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,       40.};
  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
  
  // 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;
  }
  
  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 = ptJet-rho*areaJet;
           
    Double_t varJet[6];
    varJet[0]=jet->Pt();
    varJet[1]=(Double_t)ptcorr;
    varJet[2]=(Double_t)Cent;
    varJet[3]=(Double_t)Qvec;
    varJet[4]=(Double_t)rho;
    
    AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
    if(fLeadPtMin>0 && leadTrack->Pt()<fLeadPtMin)continue;
    
    varJet[5]=(Double_t)leadTrack->Pt();
    
    ((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
  printf("AliAnalysisTaskJetSpectraAOD: Terminate() \n");
}

//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),
	88	fRejectPileup(0)
582efd60	89	{
	90	// Default constructor
	91
	92	DefineInput(0, TChain::Class());
	93	DefineOutput(1, TList::Class());
	94	DefineOutput(2, AliSpectraAODEventCuts::Class());
	95	DefineOutput(3, AliSpectraAODTrackCuts::Class());
	96
	97	}
	98
	99	//________________________________________________________________________
	100	AliAnalysisTaskJetSpectraAOD::~AliAnalysisTaskJetSpectraAOD()
	101	{
	102	delete fListJets;
	103	}
	104	//________________________________________________________________________
	105	//________________________________________________________________________
	106	void AliAnalysisTaskJetSpectraAOD::UserCreateOutputObjects()
	107	{
	108	Printf("\n\n\n\n\n\n In CreateOutput Object:");
	109
	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
	120	const Double_t ptBins[] = {0.15,5.,10.,15.,20.,25.,30.,35.,40.,50.,75.,100.,150.,200.};
	121	const Int_t nptBins=13;
	122
	123	//dimensions of THnSparse for jets
296131cf	124	const Int_t nvarjet=6;
	125	// pt_raw pt_corr cent Q vec rho pt_lead
	126	Int_t binsHistRealJet[nvarjet] = { nptBins, nptBins, fnCentBins, fnQvecBins, 40., fnptLeadBins};
	127	Double_t xminHistRealJet[nvarjet] = { 0., 0., 0., 0., 0., 0.};
	128	Double_t xmaxHistRealJet[nvarjet] = { 200., 200., 100., fQvecUpperLim, 200., 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");
	135	NSparseHistJet->SetBinEdges(1,ptBins);
	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");
	140	NSparseHistJet->GetAxis(4)->SetTitle("rho");
	141	NSparseHistJet->GetAxis(4)->SetName("rho");
296131cf	142	NSparseHistJet->GetAxis(5)->SetTitle("#it{p}_{T,lead}");
296131cf	143	NSparseHistJet->GetAxis(5)->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
	149	Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, 40.};
	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	}
	264
	265	//ESE event cuts
	266	if(!fEventCuts->IsSelected(fAOD,fTrackCuts)){
	267	if (fDebug) Printf("%s:%d ESE Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
	268	fHistEvtSelection->Fill(6);
	269	PostData(1,fOutput);
	270	PostData(2, fEventCuts);
	271	PostData(3, fTrackCuts);
	272	return;
273	}
274
275	if (fDebug) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
276	fHistEvtSelection->Fill(0.);
277	// accepted events
278	// -- end event selection --
279
582efd60	280
	281	Double_t Qvec=0.;//in case of MC we save space in the memory
	282	if(!fIsMC){
	283	if(fIsQvecCalibMode){
	284	if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
	285	else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
	286	}
	287	else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
	288	}
	289
	290	if(fIsQvecCut && (Qvec<fQvecMin \|\| Qvec>fQvecMax) ) return;
	291
	292	Double_t Cent=fEventCuts->GetCent();
296131cf	293
582efd60	294	// get background
	295	AliAODJetEventBackground* externalBackground = 0;
	296	if(fAODJets && !externalBackground && fBackgroundBranch.Length()){
	297	externalBackground = (AliAODJetEventBackground*)(fAODJets->FindListObject(fBackgroundBranch.Data()));
	298	if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
	299	}
	300
	301	Float_t rho = 0;
	302	if(externalBackground)rho = externalBackground->GetBackground(0); //default schema
	303
	304	// fetch jets
	305	TClonesArray aodJets = dynamic_cast<TClonesArray>(fAODJets->FindListObject(fJetBranchName.Data()));
	306	if(!aodJets){
	307	AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName.Data()));
	308	if(fAOD){
	309	Printf("Input AOD >>>>");
	310	fAOD->Print();
	311	}
	312	return;
	313	}
	314
	315	fListJets->Clear();
	316	for (Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
	317	AliAODJet jet = dynamic_cast<AliAODJet>((*aodJets)[iJet]);
	318	if (jet) fListJets->Add(jet);
	319	}
	320
	321
	322
	323	for(Int_t i=0; i<fListJets->GetEntries(); ++i){
	324	AliAODJet* jet = (AliAODJet*)(fListJets->At(i));
	325
	326	if((jet->Eta()<fJetEtaMin)\|\|(jet->Eta()>fJetEtaMax)) continue;
	327
	328	Double_t ptJet = jet->Pt();
	329
	330	Double_t areaJet = jet->EffectiveAreaCharged();
	331
	332	Double_t ptcorr = ptJet-rho*areaJet;
	333
296131cf	334	Double_t varJet[6];
582efd60	335	varJet[0]=jet->Pt();
	336	varJet[1]=(Double_t)ptcorr;
	337	varJet[2]=(Double_t)Cent;
	338	varJet[3]=(Double_t)Qvec;
	339	varJet[4]=(Double_t)rho;
	340
296131cf	341	AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
	342	if(fLeadPtMin>0 && leadTrack->Pt()<fLeadPtMin)continue;
	343
	344	varJet[5]=(Double_t)leadTrack->Pt();
	345
582efd60	346	((THnSparseF*)fOutput->FindObject("NSparseHistJet"))->Fill(varJet);//jet loop
	347	}
	348
	349
	350	// //track loop
	351	// for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) { //FIXME loop on aod track... should be on jet tracks???
	352	// AliAODTrack* track = fAOD->GetTrack(iTracks);
	353	// // if(!(track->TestFilterBit(1024)))continue;
	354	// if (!fTrackCuts->IsSelected(track,kTRUE)) continue;
	355	//
	356	// TH1F* h=(TH1F*)fOutput->FindObject("fHistTest");h->Fill(track->Pt());
	357	//
	358	// } // end loop on tracks
	359
	360	Double_t varEv[3];
	361	varEv[0]=Cent;
	362	varEv[1]=Qvec;
	363	varEv[2]=rho;
	364	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
	365
	366	PostData(1,fOutput);
	367	PostData(2, fEventCuts);
	368	PostData(3, fTrackCuts);
	369	//Printf("............. end of Exec");
	370
	371	}
296131cf	372
	373	//_________________________________________________________________
	374	AliVParticle AliAnalysisTaskJetSpectraAOD::LeadingTrackFromJetRefs(AliAODJet jet){
	375	if(!jet)return 0;
	376	TRefArray *refs = jet->GetRefTracks();
	377	if(!refs) return 0;
	378	AliVParticle *leading = 0;
	379	Float_t fMaxPt = 0;
	380	for(int i = 0;i<refs->GetEntriesFast();i++){
	381	AliVParticle tmp = dynamic_cast<AliVParticle>(refs->At(i));
	382	if(!tmp)continue;
	383	if(tmp->Pt()>fMaxPt){
	384	leading = tmp;
	385	fMaxPt = tmp->Pt();
	386	}
	387	}
	388	return leading;
	389	}
	390
582efd60	391	//_________________________________________________________________
	392	void AliAnalysisTaskJetSpectraAOD::Terminate(Option_t *)
	393	{
	394	// Terminate
	395	printf("AliAnalysisTaskJetSpectraAOD: Terminate() \n");
	396	}
	397
	398	//jet
	399	void AliAnalysisTaskJetSpectraAOD::SetBranchNames(const TString &branch)
	400	{
	401	fJetBranchName = branch;
	402	}
	403
	404	void AliAnalysisTaskJetSpectraAOD::SetRecBackgroundBranch(const TString &bckbranch)
	405	{
	406	fBackgroundBranch = bckbranch;
296131cf	407	}