[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetServices.cxx

// **************************************
// Task used for the correction of determiantion of reconstructed jet spectra
// Compares input (gen) and output (rec) jets   
// *******************************************


/**************************************************************************
 * 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.                  *
 **************************************************************************/

 
#include <TROOT.h>
#include <TRandom.h>
#include <TString.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
#include <TKey.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TProfile.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
#include "TDatabasePDG.h"

#include "AliAnalysisTaskJetServices.h"
#include "AliAnalysisDataContainer.h"
#include "AliAnalysisDataSlot.h"
#include "AliAnalysisManager.h"
#include "AliJetFinder.h"
#include "AliJetHeader.h"
#include "AliJetReader.h"
#include "AliJetReaderHeader.h"
#include "AliUA1JetHeaderV1.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliAODTrack.h"
#include "AliAODJet.h"
#include "AliAODMCParticle.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliGenPythiaEventHeader.h"
#include "AliJetKineReaderHeader.h"
#include "AliGenCocktailEventHeader.h"
#include "AliInputEventHandler.h"
#include "AliPhysicsSelection.h"


#include "AliAnalysisHelperJetTasks.h"

ClassImp(AliAnalysisTaskJetServices)

AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
  fUseAODInput(kFALSE),
  fUsePhysicsSelection(kFALSE),
  fRealData(kFALSE),
  fSelectionInfoESD(0),
  fAvgTrials(1),
  fVtxXMean(0),
  fVtxYMean(0),
  fVtxZMean(0),
  fVtxRCut(1.),
  fVtxZCut(8.),
  fh1Xsec(0x0),
  fh1Trials(0x0),
  fh1PtHard(0x0),
  fh1PtHardTrials(0x0),
  fh1SelectionInfoESD(0x0),
  fh2TriggerCount(0x0),
  fh2ESDTriggerCount(0x0),
  fh2TriggerVtx(0x0),
  fh2ESDTriggerVtx(0x0),
  fh2ESDTriggerRun(0x0),
  fh2VtxXY(0x0),
  fHistList(0x0)  
{
  fRunRange[0] = fRunRange[1] = 0; 
}

AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
  AliAnalysisTaskSE(name),
  fUseAODInput(kFALSE),
  fUsePhysicsSelection(kFALSE),
  fRealData(kFALSE),
  fSelectionInfoESD(0),
  fAvgTrials(1),
  fVtxXMean(0),
  fVtxYMean(0),
  fVtxZMean(0),
  fVtxRCut(1.),
  fVtxZCut(8.),
  fh1Xsec(0x0),
  fh1Trials(0x0),
  fh1PtHard(0x0),
  fh1PtHardTrials(0x0),
  fh1SelectionInfoESD(0x0),
  fh2TriggerCount(0x0),
  fh2ESDTriggerCount(0x0),
  fh2TriggerVtx(0x0),
  fh2ESDTriggerVtx(0x0),
  fh2ESDTriggerRun(0x0),
  fh2VtxXY(0x0),
  fHistList(0x0)  
{
  fRunRange[0] = fRunRange[1] = 0; 
  DefineOutput(1,TList::Class());
}


Bool_t AliAnalysisTaskJetServices::Notify()
{
  //
  // Implemented Notify() to read the cross sections
  // and number of trials from pyxsec.root
  // 

  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
  Float_t xsection = 0;
  Float_t ftrials  = 1;

  fAvgTrials = 1;
  if(tree){
    TFile *curfile = tree->GetCurrentFile();
    if (!curfile) {
      Error("Notify","No current file");
      return kFALSE;
    }
    if(!fh1Xsec||!fh1Trials){
      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
      return kFALSE;
    }
    AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
    fh1Xsec->Fill("<#sigma>",xsection);
    // construct a poor man average trials 
    Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
    if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
  }  
  return kTRUE;
}

void AliAnalysisTaskJetServices::UserCreateOutputObjects()
{

  //
  // Create the output container
  //


  if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");

  OpenFile(1);
  if(!fHistList)fHistList = new TList();

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);
  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
  fHistList->Add(fh1Xsec);

  fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
  fHistList->Add(fh1Trials);

  const Int_t nBinPt = 100;
  Double_t binLimitsPt[nBinPt+1];
  for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
    if(iPt == 0){
      binLimitsPt[iPt] = 0.0;
    }
    else {// 1.0
      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
    }
  }
  
  fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
  fHistList->Add(fh2TriggerCount);

  fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
  fHistList->Add(fh2ESDTriggerCount);

  fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
  fHistList->Add(fh2TriggerVtx);

  fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
  fHistList->Add(fh2ESDTriggerVtx);
  

  fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
  fHistList->Add(fh1PtHard);
  fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
  fHistList->Add(fh1PtHardTrials);
  fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
				 AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);

  fHistList->Add(fh1SelectionInfoESD);
  // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range  
  // 3 triggers BB BE/EB EE

  fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Trigger vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
  fHistList->Add(fh2ESDTriggerRun);

  fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
  fHistList->Add(fh2VtxXY);
  // =========== Switch on Sumw2 for all histos ===========
  for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
    TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
    if (h1){
      h1->Sumw2();
      continue;
    }
    THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
    if(hn)hn->Sumw2();
  }


  TH1::AddDirectory(oldStatus);
}

void AliAnalysisTaskJetServices::Init()
{
  //
  // Initialization
  //
  if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");

}

void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
{

  //
  // Execute analysis for current event
  //
 
  AliAODEvent *aod = 0;
  AliESDEvent *esd = 0;

  AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
  fSelectionInfoESD = 0; // reset
  AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false


  if(fUseAODInput){    
    aod = dynamic_cast<AliAODEvent*>(InputEvent());
    if(!aod){
      Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
      return;
    }
    // fethc the header
  }
  else{
    //  assume that the AOD is in the general output...
    aod  = AODEvent();
    if(!aod){
      Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
      return;
    }
    esd = dynamic_cast<AliESDEvent*>(InputEvent());
  }
  
  fSelectionInfoESD |= AliAnalysisHelperJetTasks::kNone;

  if(esd){
    Float_t run = (Float_t)esd->GetRunNumber();
    const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
    Float_t zvtx = -999;
    Float_t xvtx = -999;
    Float_t yvtx = -999;

    if(vtxESD->GetNContributors()>2){
      zvtx = vtxESD->GetZ();
      yvtx = vtxESD->GetY();
      xvtx = vtxESD->GetX();
    }

    Int_t iTrig = -1;
    if(esd->GetFiredTriggerClasses().Contains("CINT1B")
       ||esd->GetFiredTriggerClasses().Contains("CSMBB")
       ||esd->GetFiredTriggerClasses().Contains("MB1")
       ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
      iTrig = 0;
      fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchBunch;
    }
    else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
	    ||esd->GetFiredTriggerClasses().Contains("CSMBA")
	    ||esd->GetFiredTriggerClasses().Contains("CINT6A")
	    ||esd->GetFiredTriggerClasses().Contains("CINT1C")
	    ||esd->GetFiredTriggerClasses().Contains("CSMBC")
	    ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
      // empty bunch or bunch empty
      iTrig = 1;
      fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchEmpty;
    }
    else if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
       ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
      iTrig = 2;
      fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kEmptyEmpty;
    }

    
    if(iTrig>=0){
      iTrig *= 3;
      fh2ESDTriggerRun->Fill(run,iTrig+1);
      if(vtxESD->GetNContributors()>2){
	fh2ESDTriggerRun->Fill(run,iTrig+2);
	fh2VtxXY->Fill(xvtx,yvtx);
      }
      xvtx -= fVtxXMean; 
      yvtx -= fVtxYMean; 
      zvtx -= fVtxZMean; 
      Float_t r2 = xvtx *xvtx + yvtx *yvtx; 
      if(TMath::Abs(zvtx)<fVtxZCut&&r2<fVtxRCut)fh2ESDTriggerRun->Fill(run,iTrig+3);
    }
    else{
      fh2ESDTriggerRun->Fill(run,0);
    }


  }
  

  // Apply additional constraints
  Bool_t esdEventSelected = IsEventSelectedESD(esd);
  Bool_t esdEventPileUp = IsEventPileUpESD(esd);
  Bool_t esdEventCosmic = IsEventCosmicESD(esd);

  Bool_t aodEventSelected = IsEventSelectedAOD(aod);

  Bool_t physicsSelection = fInputHandler->IsEventSelected();

  if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
  if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
  if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
  if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;


  // here we have all selection information, fill histogram
  for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
    if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
  }
  AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false

  if(esd&&aod&&fDebug){
    if(esdEventSelected&&!aodEventSelected){
      Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
      const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
      const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
      Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
      vtxESD->Print();
      Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
      vtxAOD->Print();
    }
  }

  // loop over all possible triggers for esd 

  if(esd){
    const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
    //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
    TString vtxName(vtxESD->GetName());
    for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
      Bool_t esdTrig = kFALSE;
      esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
      if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
      Bool_t cand = fInputHandler->IsEventSelected();
      if(cand){
	fh2ESDTriggerCount->Fill(it,kSelectedALICE); 
      }
      if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
      if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){
	if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex);
	Float_t zvtx = vtxESD->GetZ();
	if(esdEventSelected&&esdTrig){
	  fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
	  fh2ESDTriggerVtx->Fill(it,zvtx);
	}
      }
      if(cand&&esdEventSelected){
	fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
	fh2ESDTriggerCount->Fill(it,kSelected);
	AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
      }
    }
  }

  if(aod){
    const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
    //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
    TString vtxTitle(vtxAOD->GetTitle());
    for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
      Bool_t aodTrig = kFALSE;
      aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
      if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
      if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
	if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
	Float_t zvtx = vtxAOD->GetZ();
	if(aodTrig&&aodEventSelected){
	  fh2TriggerVtx->Fill(it,zvtx);
	  fh2TriggerCount->Fill(it,kTriggeredVertexIn);
	}
      }
    }
  }

  if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);

  
  Double_t ptHard = 0; 
  Double_t nTrials = 1; // Trials for MC trigger 

  fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
  AliMCEvent* mcEvent = MCEvent();
  //    AliStack *pStack = 0; 
  if(mcEvent){
    AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
    if(pythiaGenHeader){
      nTrials = pythiaGenHeader->Trials();
      ptHard  = pythiaGenHeader->GetPtHard();
      int iProcessType = pythiaGenHeader->ProcessType();
      // 11 f+f -> f+f
      // 12 f+barf -> f+barf
      // 13 f+barf -> g+g
      // 28 f+g -> f+g
      // 53 g+g -> f+barf
      // 68 g+g -> g+g
      if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
      if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
      fh1PtHard->Fill(ptHard);
      fh1PtHardTrials->Fill(ptHard,nTrials);

    }// if pythia gen header
  }

  // trigger selection
  

  PostData(1, fHistList);
}

Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
  if(!esd)return kFALSE;
  const AliESDVertex *vtx = esd->GetPrimaryVertex();
  // We can check the type of the vertex by its name and title
  // if vertexer failed title is empty (default c'tor) and ncontributors is 0
  // vtx       name                  title
  // Tracks    PrimaryVertex         VertexerTracksNoConstraint
  // TPC       TPCVertex             VertexerTracksNoConstraint
  // SPD       SPDVertex             vertexer: 3D
  // SPD       SPDVertex             vertexer: Z
  

  if(vtx->GetNContributors()<3)return kFALSE;

  // do not want tpc only primary vertex
  TString vtxName(vtx->GetName());
  if(vtxName.Contains("TPCVertex"))return kFALSE;
  Float_t zvtx = vtx->GetZ();
  Float_t yvtx = vtx->GetY();
  Float_t xvtx = vtx->GetX();

  // here we may fill the histos for run dependence....

  xvtx -= fVtxXMean;
  yvtx -= fVtxYMean;
  zvtx -= fVtxZMean;

  Float_t r2   = yvtx*yvtx+xvtx*xvtx;      

  Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
  return eventSel;
}

Bool_t AliAnalysisTaskJetServices::IsEventPileUpESD(AliESDEvent* esd){
  if(!esd)return kFALSE;
  return esd->IsPileupFromSPD();
}

Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
  if(!esd)return kFALSE;
  // add track cuts for which we look for cosmics...
  return kTRUE;
}


Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
  if(!aod)return kFALSE;
  const AliAODVertex *vtx = aod->GetPrimaryVertex();
  // We can check the type of the vertex by its name and title
  // if vertexer failed title is empty (default c'tor) and ncontributors is 0
  // vtx       name                  title
  // Tracks    PrimaryVertex         VertexerTracksNoConstraint
  // TPC       TPCVertex             VertexerTracksNoConstraint
  // SPD       SPDVertex             vertexer: 3D
  // SPD       SPDVertex             vertexer: Z
  

  if(vtx->GetNContributors()<3)return kFALSE;

  // do not want tpc only primary vertex
  TString vtxName(vtx->GetName());
  if(vtxName.Contains("TPCVertex"))return kFALSE;

  Float_t zvtx = vtx->GetZ();
  Float_t yvtx = vtx->GetY();
  Float_t xvtx = vtx->GetX();

  // here we may fill the histos for run dependence....

  xvtx -= fVtxXMean;
  yvtx -= fVtxYMean;
  zvtx -= fVtxZMean;

  Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
  Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
  return eventSel;
}


void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
}
Commit	Line	Data
bf7b8731	1	// **************************************
	2	// Task used for the correction of determiantion of reconstructed jet spectra
	3	// Compares input (gen) and output (rec) jets
	4	// *******************************************
	5
	6
	7	/**************************************************************************
	8	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	9	* *
	10	* Author: The ALICE Off-line Project. *
	11	* Contributors are mentioned in the code where appropriate. *
	12	* *
	13	* Permission to use, copy, modify and distribute this software and its *
	14	* documentation strictly for non-commercial purposes is hereby granted *
	15	* without fee, provided that the above copyright notice appears in all *
	16	* copies and that both the copyright notice and this permission notice *
	17	* appear in the supporting documentation. The authors make no claims *
	18	* about the suitability of this software for any purpose. It is *
	19	* provided "as is" without express or implied warranty. *
	20	**************************************************************************/
	21
	22
	23	#include <TROOT.h>
	24	#include <TRandom.h>
b1cd90b5	25	#include <TString.h>
bf7b8731	26	#include <TSystem.h>
	27	#include <TInterpreter.h>
	28	#include <TChain.h>
	29	#include <TFile.h>
	30	#include <TKey.h>
	31	#include <TH1F.h>
	32	#include <TH2F.h>
	33	#include <TH3F.h>
	34	#include <TProfile.h>
	35	#include <TList.h>
	36	#include <TLorentzVector.h>
	37	#include <TClonesArray.h>
c0997643	38	#include "TDatabasePDG.h"
bf7b8731	39
bf7b8731	40	#include "AliAnalysisTaskJetServices.h"
fe669ac6	41	#include "AliAnalysisDataContainer.h"
fe669ac6	42	#include "AliAnalysisDataSlot.h"
bf7b8731	43	#include "AliAnalysisManager.h"
	44	#include "AliJetFinder.h"
	45	#include "AliJetHeader.h"
	46	#include "AliJetReader.h"
	47	#include "AliJetReaderHeader.h"
	48	#include "AliUA1JetHeaderV1.h"
	49	#include "AliESDEvent.h"
	50	#include "AliAODEvent.h"
	51	#include "AliAODHandler.h"
	52	#include "AliAODTrack.h"
	53	#include "AliAODJet.h"
	54	#include "AliAODMCParticle.h"
	55	#include "AliMCEventHandler.h"
	56	#include "AliMCEvent.h"
	57	#include "AliStack.h"
	58	#include "AliGenPythiaEventHeader.h"
	59	#include "AliJetKineReaderHeader.h"
	60	#include "AliGenCocktailEventHeader.h"
	61	#include "AliInputEventHandler.h"
fe669ac6	62	#include "AliPhysicsSelection.h"
bf7b8731	63
	64
	65	#include "AliAnalysisHelperJetTasks.h"
	66
	67	ClassImp(AliAnalysisTaskJetServices)
	68
	69	AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
	70	fUseAODInput(kFALSE),
cc0649e4	71	fUsePhysicsSelection(kFALSE),
cc0649e4	72	fRealData(kFALSE),
f2dd0695	73	fSelectionInfoESD(0),
bf7b8731	74	fAvgTrials(1),
c0997643	75	fVtxXMean(0),
	76	fVtxYMean(0),
	77	fVtxZMean(0),
	78	fVtxRCut(1.),
	79	fVtxZCut(8.),
bf7b8731	80	fh1Xsec(0x0),
	81	fh1Trials(0x0),
	82	fh1PtHard(0x0),
	83	fh1PtHardTrials(0x0),
f2dd0695	84	fh1SelectionInfoESD(0x0),
bf7b8731	85	fh2TriggerCount(0x0),
	86	fh2ESDTriggerCount(0x0),
	87	fh2TriggerVtx(0x0),
	88	fh2ESDTriggerVtx(0x0),
b5a3f310	89	fh2ESDTriggerRun(0x0),
b5a3f310	90	fh2VtxXY(0x0),
bf7b8731	91	fHistList(0x0)
bf7b8731	92	{
b5a3f310	93	fRunRange[0] = fRunRange[1] = 0;
bf7b8731	94	}
	95
	96	AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
	97	AliAnalysisTaskSE(name),
	98	fUseAODInput(kFALSE),
cc0649e4	99	fUsePhysicsSelection(kFALSE),
cc0649e4	100	fRealData(kFALSE),
f2dd0695	101	fSelectionInfoESD(0),
bf7b8731	102	fAvgTrials(1),
c0997643	103	fVtxXMean(0),
	104	fVtxYMean(0),
	105	fVtxZMean(0),
	106	fVtxRCut(1.),
	107	fVtxZCut(8.),
bf7b8731	108	fh1Xsec(0x0),
	109	fh1Trials(0x0),
	110	fh1PtHard(0x0),
	111	fh1PtHardTrials(0x0),
f2dd0695	112	fh1SelectionInfoESD(0x0),
bf7b8731	113	fh2TriggerCount(0x0),
	114	fh2ESDTriggerCount(0x0),
	115	fh2TriggerVtx(0x0),
	116	fh2ESDTriggerVtx(0x0),
b5a3f310	117	fh2ESDTriggerRun(0x0),
b5a3f310	118	fh2VtxXY(0x0),
bf7b8731	119	fHistList(0x0)
bf7b8731	120	{
b5a3f310	121	fRunRange[0] = fRunRange[1] = 0;
bf7b8731	122	DefineOutput(1,TList::Class());
	123	}
	124
	125
	126
	127	Bool_t AliAnalysisTaskJetServices::Notify()
	128	{
	129	//
	130	// Implemented Notify() to read the cross sections
	131	// and number of trials from pyxsec.root
	132	//
	133
	134	TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
	135	Float_t xsection = 0;
	136	Float_t ftrials = 1;
	137
	138	fAvgTrials = 1;
	139	if(tree){
	140	TFile *curfile = tree->GetCurrentFile();
	141	if (!curfile) {
	142	Error("Notify","No current file");
	143	return kFALSE;
	144	}
	145	if(!fh1Xsec\|\|!fh1Trials){
	146	Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
	147	return kFALSE;
	148	}
	149	AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
	150	fh1Xsec->Fill("<#sigma>",xsection);
	151	// construct a poor man average trials
	152	Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
	153	if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
	154	}
	155	return kTRUE;
	156	}
	157
	158	void AliAnalysisTaskJetServices::UserCreateOutputObjects()
	159	{
	160
	161	//
	162	// Create the output container
	163	//
	164
	165
bf7b8731	166	if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
	167
	168	OpenFile(1);
	169	if(!fHistList)fHistList = new TList();
	170
	171	Bool_t oldStatus = TH1::AddDirectoryStatus();
	172	TH1::AddDirectory(kFALSE);
bf7b8731	173	fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
	174	fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
	175	fHistList->Add(fh1Xsec);
	176
	177	fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
	178	fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
	179	fHistList->Add(fh1Trials);
	180
	181	const Int_t nBinPt = 100;
	182	Double_t binLimitsPt[nBinPt+1];
	183	for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
	184	if(iPt == 0){
	185	binLimitsPt[iPt] = 0.0;
	186	}
	187	else {// 1.0
	188	binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
	189	}
	190	}
	191
	192	fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
	193	fHistList->Add(fh2TriggerCount);
	194
	195	fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
	196	fHistList->Add(fh2ESDTriggerCount);
	197
	198	fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
	199	fHistList->Add(fh2TriggerVtx);
	200
	201	fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
	202	fHistList->Add(fh2ESDTriggerVtx);
	203
	204
	205	fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
	206	fHistList->Add(fh1PtHard);
	207	fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
	208	fHistList->Add(fh1PtHardTrials);
f2dd0695	209	fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
f2dd0695	210	AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
bf7b8731	211
f2dd0695	212	fHistList->Add(fh1SelectionInfoESD);
b5a3f310	213	// 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
	214	// 3 triggers BB BE/EB EE
	215
	216	fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Trigger vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
	217	fHistList->Add(fh2ESDTriggerRun);
	218
	219	fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
	220	fHistList->Add(fh2VtxXY);
bf7b8731	221	// =========== Switch on Sumw2 for all histos ===========
	222	for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
	223	TH1 h1 = dynamic_cast<TH1>(fHistList->At(i));
	224	if (h1){
	225	h1->Sumw2();
	226	continue;
	227	}
	228	THnSparse hn = dynamic_cast<THnSparse>(fHistList->At(i));
	229	if(hn)hn->Sumw2();
	230	}
	231
	232
	233	TH1::AddDirectory(oldStatus);
	234	}
	235
	236	void AliAnalysisTaskJetServices::Init()
	237	{
	238	//
	239	// Initialization
	240	//
bf7b8731	241	if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
	242
	243	}
	244
	245	void AliAnalysisTaskJetServices::UserExec(Option_t /option*/)
	246	{
	247
	248	//
	249	// Execute analysis for current event
	250	//
	251
	252	AliAODEvent *aod = 0;
	253	AliESDEvent *esd = 0;
	254
	255	AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
f2dd0695	256	fSelectionInfoESD = 0; // reset
	257	AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
	258
bf7b8731	259
	260	if(fUseAODInput){
	261	aod = dynamic_cast<AliAODEvent*>(InputEvent());
	262	if(!aod){
	263	Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
	264	return;
	265	}
	266	// fethc the header
	267	}
	268	else{
	269	// assume that the AOD is in the general output...
	270	aod = AODEvent();
	271	if(!aod){
	272	Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
	273	return;
	274	}
	275	esd = dynamic_cast<AliESDEvent*>(InputEvent());
	276	}
	277
f2dd0695	278	fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kNone;
f2dd0695	279
b5a3f310	280	if(esd){
	281	Float_t run = (Float_t)esd->GetRunNumber();
	282	const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
	283	Float_t zvtx = -999;
	284	Float_t xvtx = -999;
	285	Float_t yvtx = -999;
	286
3568c3d6	287	if(vtxESD->GetNContributors()>2){
b5a3f310	288	zvtx = vtxESD->GetZ();
	289	yvtx = vtxESD->GetY();
	290	xvtx = vtxESD->GetX();
	291	}
	292
	293	Int_t iTrig = -1;
	294	if(esd->GetFiredTriggerClasses().Contains("CINT1B")
	295	\|\|esd->GetFiredTriggerClasses().Contains("CSMBB")
	296	\|\|esd->GetFiredTriggerClasses().Contains("MB1")
	297	\|\|esd->GetFiredTriggerClasses().Contains("CINT6B")){
	298	iTrig = 0;
f2dd0695	299	fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kBunchBunch;
b5a3f310	300	}
	301	else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
	302	\|\|esd->GetFiredTriggerClasses().Contains("CSMBA")
	303	\|\|esd->GetFiredTriggerClasses().Contains("CINT6A")
	304	\|\|esd->GetFiredTriggerClasses().Contains("CINT1C")
	305	\|\|esd->GetFiredTriggerClasses().Contains("CSMBC")
	306	\|\|esd->GetFiredTriggerClasses().Contains("CINT6C")){
	307	// empty bunch or bunch empty
	308	iTrig = 1;
f2dd0695	309	fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kBunchEmpty;
b5a3f310	310	}
f2dd0695	311	else if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
b5a3f310	312	\|\|esd->GetFiredTriggerClasses().Contains("CINT6-E")){
b5a3f310	313	iTrig = 2;
f2dd0695	314	fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kEmptyEmpty;
b5a3f310	315	}
	316
	317
	318	if(iTrig>=0){
	319	iTrig *= 3;
	320	fh2ESDTriggerRun->Fill(run,iTrig+1);
3568c3d6	321	if(vtxESD->GetNContributors()>2){
b5a3f310	322	fh2ESDTriggerRun->Fill(run,iTrig+2);
	323	fh2VtxXY->Fill(xvtx,yvtx);
	324	}
c0997643	325	xvtx -= fVtxXMean;
	326	yvtx -= fVtxYMean;
	327	zvtx -= fVtxZMean;
cf049e94	328	Float_t r2 = xvtx xvtx + yvtx yvtx;
c0997643	329	if(TMath::Abs(zvtx)<fVtxZCut&&r2<fVtxRCut)fh2ESDTriggerRun->Fill(run,iTrig+3);
b5a3f310	330	}
	331	else{
	332	fh2ESDTriggerRun->Fill(run,0);
	333	}
	334
	335
	336	}
bf7b8731	337
bf7b8731	338
c0997643	339	// Apply additional constraints
c0997643	340	Bool_t esdEventSelected = IsEventSelectedESD(esd);
f2dd0695	341	Bool_t esdEventPileUp = IsEventPileUpESD(esd);
	342	Bool_t esdEventCosmic = IsEventCosmicESD(esd);
	343
c0997643	344	Bool_t aodEventSelected = IsEventSelectedAOD(aod);
f2dd0695	345
	346	Bool_t physicsSelection = fInputHandler->IsEventSelected();
	347
	348	if(esdEventSelected) fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kVertexIn;
	349	if(esdEventPileUp) fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kIsPileUp;
	350	if(esdEventCosmic) fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kIsCosmic;
	351	if(physicsSelection) fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kPhysicsSelection;
	352
	353
	354	// here we have all selection information, fill histogram
	355	for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
	356	if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
	357	}
	358	AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
c0997643	359
	360	if(esd&&aod&&fDebug){
	361	if(esdEventSelected&&!aodEventSelected){
	362	Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
	363	const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
bf7b8731	364	const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
c0997643	365	Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
	366	vtxESD->Print();
	367	Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
	368	vtxAOD->Print();
bf7b8731	369	}
c0997643	370	}
	371
	372	// loop over all possible triggers for esd
	373
	374	if(esd){
	375	const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
	376	// Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
	377	TString vtxName(vtxESD->GetName());
	378	for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
	379	Bool_t esdTrig = kFALSE;
	380	esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
	381	if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
216601f0	382	Bool_t cand = fInputHandler->IsEventSelected();
	383	if(cand){
	384	fh2ESDTriggerCount->Fill(it,kSelectedALICE);
216601f0	385	}
cc0649e4	386	if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
3568c3d6	387	if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){
c0997643	388	if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex);
bf7b8731	389	Float_t zvtx = vtxESD->GetZ();
c0997643	390	if(esdEventSelected&&esdTrig){
bf7b8731	391	fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
cc89bb69	392	fh2ESDTriggerVtx->Fill(it,zvtx);
bf7b8731	393	}
bf7b8731	394	}
c0997643	395	if(cand&&esdEventSelected){
	396	fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
	397	fh2ESDTriggerCount->Fill(it,kSelected);
	398	AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
	399	}
bf7b8731	400	}
bf7b8731	401	}
bf7b8731	402
c0997643	403	if(aod){
	404	const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
	405	// Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
	406	TString vtxTitle(vtxAOD->GetTitle());
	407	for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
	408	Bool_t aodTrig = kFALSE;
	409	aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
	410	if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
3568c3d6	411	if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
c0997643	412	if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
c0997643	413	Float_t zvtx = vtxAOD->GetZ();
c0997643	414	if(aodTrig&&aodEventSelected){
cc89bb69	415	fh2TriggerVtx->Fill(it,zvtx);
c0997643	416	fh2TriggerCount->Fill(it,kTriggeredVertexIn);
	417	}
	418	}
	419	}
	420	}
bf7b8731	421
	422	if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
	423
	424
	425	Double_t ptHard = 0;
	426	Double_t nTrials = 1; // Trials for MC trigger
	427
	428	fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
	429	AliMCEvent* mcEvent = MCEvent();
	430	// AliStack *pStack = 0;
	431	if(mcEvent){
	432	AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
	433	if(pythiaGenHeader){
	434	nTrials = pythiaGenHeader->Trials();
	435	ptHard = pythiaGenHeader->GetPtHard();
	436	int iProcessType = pythiaGenHeader->ProcessType();
	437	// 11 f+f -> f+f
	438	// 12 f+barf -> f+barf
	439	// 13 f+barf -> g+g
	440	// 28 f+g -> f+g
	441	// 53 g+g -> f+barf
	442	// 68 g+g -> g+g
	443	if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
	444	if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
	445	fh1PtHard->Fill(ptHard);
	446	fh1PtHardTrials->Fill(ptHard,nTrials);
	447
	448	}// if pythia gen header
	449	}
	450
	451	// trigger selection
	452
	453
	454	PostData(1, fHistList);
	455	}
	456
c0997643	457	Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
	458	if(!esd)return kFALSE;
	459	const AliESDVertex *vtx = esd->GetPrimaryVertex();
	460	// We can check the type of the vertex by its name and title
	461	// if vertexer failed title is empty (default c'tor) and ncontributors is 0
	462	// vtx name title
	463	// Tracks PrimaryVertex VertexerTracksNoConstraint
	464	// TPC TPCVertex VertexerTracksNoConstraint
	465	// SPD SPDVertex vertexer: 3D
	466	// SPD SPDVertex vertexer: Z
	467
	468
	469
3568c3d6	470	if(vtx->GetNContributors()<3)return kFALSE;
c0997643	471
	472	// do not want tpc only primary vertex
	473	TString vtxName(vtx->GetName());
	474	if(vtxName.Contains("TPCVertex"))return kFALSE;
	475	Float_t zvtx = vtx->GetZ();
	476	Float_t yvtx = vtx->GetY();
	477	Float_t xvtx = vtx->GetX();
	478
	479	// here we may fill the histos for run dependence....
	480
	481	xvtx -= fVtxXMean;
	482	yvtx -= fVtxYMean;
	483	zvtx -= fVtxZMean;
	484
	485	Float_t r2 = yvtxyvtx+xvtxxvtx;
	486
	487	Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
	488	return eventSel;
	489	}
	490
f2dd0695	491	Bool_t AliAnalysisTaskJetServices::IsEventPileUpESD(AliESDEvent* esd){
	492	if(!esd)return kFALSE;
	493	return esd->IsPileupFromSPD();
	494	}
	495
	496	Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
	497	if(!esd)return kFALSE;
	498	// add track cuts for which we look for cosmics...
	499	return kTRUE;
	500	}
	501
	502
c0997643	503	Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
	504	if(!aod)return kFALSE;
	505	const AliAODVertex *vtx = aod->GetPrimaryVertex();
	506	// We can check the type of the vertex by its name and title
	507	// if vertexer failed title is empty (default c'tor) and ncontributors is 0
	508	// vtx name title
	509	// Tracks PrimaryVertex VertexerTracksNoConstraint
	510	// TPC TPCVertex VertexerTracksNoConstraint
	511	// SPD SPDVertex vertexer: 3D
	512	// SPD SPDVertex vertexer: Z
	513
	514
	515
3568c3d6	516	if(vtx->GetNContributors()<3)return kFALSE;
c0997643	517
	518	// do not want tpc only primary vertex
	519	TString vtxName(vtx->GetName());
	520	if(vtxName.Contains("TPCVertex"))return kFALSE;
	521
	522	Float_t zvtx = vtx->GetZ();
	523	Float_t yvtx = vtx->GetY();
	524	Float_t xvtx = vtx->GetX();
	525
	526	// here we may fill the histos for run dependence....
	527
	528	xvtx -= fVtxXMean;
	529	yvtx -= fVtxYMean;
	530	zvtx -= fVtxZMean;
	531
	532	Float_t r2 = yvtxyvtx+xvtxxvtx;
	533	Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
	534	return eventSel;
	535	}
	536
	537
bf7b8731	538
f2dd0695	539
bf7b8731	540	void AliAnalysisTaskJetServices::Terminate(Option_t /option*/)
	541	{
	542	// Terminate analysis
	543	//
bf7b8731	544	}