[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.),
  fPtMinCosmic(5.),
  fRIsolMinCosmic(3.),
  fMaxCosmicAngle(0.01),
  fh1Xsec(0x0),
  fh1Trials(0x0),
  fh1PtHard(0x0),
  fh1PtHardTrials(0x0),
  fh1SelectionInfoESD(0x0),
  fh2TriggerCount(0x0),
  fh2ESDTriggerCount(0x0),
  fh2TriggerVtx(0x0),
  fh2ESDTriggerVtx(0x0),
  fh2ESDTriggerRun(0x0),
  fh2VtxXY(0x0),
  fh1NCosmicsPerEvent(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.),
  fPtMinCosmic(5.),
  fRIsolMinCosmic(3.),
  fMaxCosmicAngle(0.01),
  fh1Xsec(0x0),
  fh1Trials(0x0),
  fh1PtHard(0x0),
  fh1PtHardTrials(0x0),
  fh1SelectionInfoESD(0x0),
  fh2TriggerCount(0x0),
  fh2ESDTriggerCount(0x0),
  fh2TriggerVtx(0x0),
  fh2ESDTriggerVtx(0x0),
  fh2ESDTriggerRun(0x0),
  fh2VtxXY(0x0),
  fh1NCosmicsPerEvent(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();
  }

  fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
  fHistList->Add(fh1NCosmicsPerEvent),


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

  Bool_t isCosmic = kFALSE;
  Int_t nTracks = esd->GetNumberOfTracks();
  Int_t nCosmicCandidates = 0;

  for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
    AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
    if (!track1)  continue;
    if(track1->Pt()<fPtMinCosmic) continue;
    //Start 2nd track loop to look for correlations
    for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
      AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
      if(!track2) continue;
      if(track2->Pt()<fPtMinCosmic) continue;
      //Check if back-to-back
      Double_t mom1[3],mom2[3];
      track1->GetPxPyPz(mom1);
      track2->GetPxPyPz(mom2);
      TVector3 momv1(mom1[0],mom1[1],mom1[2]);
      TVector3 momv2(mom2[0],mom2[1],mom2[2]);
      Float_t theta = (float)(momv1.Phi()-momv2.Phi());
      if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();

      Float_t deltaPhi = track1->Phi()-track2->Phi();
      if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();

      Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
      if(rIsol<fRIsolMinCosmic) continue;

      if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
	nCosmicCandidates+=1.;
	isCosmic = kTRUE;
      }
      
    }
  }

  fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);

  return isCosmic;
}


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.),
df7848fc	80	fPtMinCosmic(5.),
	81	fRIsolMinCosmic(3.),
	82	fMaxCosmicAngle(0.01),
bf7b8731	83	fh1Xsec(0x0),
	84	fh1Trials(0x0),
	85	fh1PtHard(0x0),
	86	fh1PtHardTrials(0x0),
f2dd0695	87	fh1SelectionInfoESD(0x0),
bf7b8731	88	fh2TriggerCount(0x0),
	89	fh2ESDTriggerCount(0x0),
	90	fh2TriggerVtx(0x0),
	91	fh2ESDTriggerVtx(0x0),
b5a3f310	92	fh2ESDTriggerRun(0x0),
b5a3f310	93	fh2VtxXY(0x0),
df7848fc	94	fh1NCosmicsPerEvent(0x0),
bf7b8731	95	fHistList(0x0)
bf7b8731	96	{
b5a3f310	97	fRunRange[0] = fRunRange[1] = 0;
bf7b8731	98	}
	99
	100	AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
	101	AliAnalysisTaskSE(name),
	102	fUseAODInput(kFALSE),
cc0649e4	103	fUsePhysicsSelection(kFALSE),
cc0649e4	104	fRealData(kFALSE),
f2dd0695	105	fSelectionInfoESD(0),
bf7b8731	106	fAvgTrials(1),
c0997643	107	fVtxXMean(0),
	108	fVtxYMean(0),
	109	fVtxZMean(0),
	110	fVtxRCut(1.),
	111	fVtxZCut(8.),
df7848fc	112	fPtMinCosmic(5.),
	113	fRIsolMinCosmic(3.),
	114	fMaxCosmicAngle(0.01),
bf7b8731	115	fh1Xsec(0x0),
	116	fh1Trials(0x0),
	117	fh1PtHard(0x0),
	118	fh1PtHardTrials(0x0),
f2dd0695	119	fh1SelectionInfoESD(0x0),
bf7b8731	120	fh2TriggerCount(0x0),
	121	fh2ESDTriggerCount(0x0),
	122	fh2TriggerVtx(0x0),
	123	fh2ESDTriggerVtx(0x0),
b5a3f310	124	fh2ESDTriggerRun(0x0),
b5a3f310	125	fh2VtxXY(0x0),
df7848fc	126	fh1NCosmicsPerEvent(0x0),
bf7b8731	127	fHistList(0x0)
bf7b8731	128	{
b5a3f310	129	fRunRange[0] = fRunRange[1] = 0;
bf7b8731	130	DefineOutput(1,TList::Class());
	131	}
	132
	133
	134
	135	Bool_t AliAnalysisTaskJetServices::Notify()
	136	{
	137	//
	138	// Implemented Notify() to read the cross sections
	139	// and number of trials from pyxsec.root
	140	//
	141
	142	TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
	143	Float_t xsection = 0;
	144	Float_t ftrials = 1;
	145
	146	fAvgTrials = 1;
	147	if(tree){
	148	TFile *curfile = tree->GetCurrentFile();
	149	if (!curfile) {
	150	Error("Notify","No current file");
	151	return kFALSE;
	152	}
	153	if(!fh1Xsec\|\|!fh1Trials){
	154	Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
	155	return kFALSE;
	156	}
	157	AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
	158	fh1Xsec->Fill("<#sigma>",xsection);
	159	// construct a poor man average trials
	160	Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
	161	if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
	162	}
	163	return kTRUE;
	164	}
	165
	166	void AliAnalysisTaskJetServices::UserCreateOutputObjects()
	167	{
	168
	169	//
	170	// Create the output container
	171	//
	172
	173
bf7b8731	174	if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
	175
	176	OpenFile(1);
	177	if(!fHistList)fHistList = new TList();
	178
	179	Bool_t oldStatus = TH1::AddDirectoryStatus();
	180	TH1::AddDirectory(kFALSE);
bf7b8731	181	fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
	182	fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
	183	fHistList->Add(fh1Xsec);
	184
	185	fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
	186	fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
	187	fHistList->Add(fh1Trials);
	188
	189	const Int_t nBinPt = 100;
	190	Double_t binLimitsPt[nBinPt+1];
	191	for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
	192	if(iPt == 0){
	193	binLimitsPt[iPt] = 0.0;
	194	}
	195	else {// 1.0
	196	binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
	197	}
	198	}
	199
	200	fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
	201	fHistList->Add(fh2TriggerCount);
	202
	203	fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
	204	fHistList->Add(fh2ESDTriggerCount);
	205
	206	fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
	207	fHistList->Add(fh2TriggerVtx);
	208
	209	fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
	210	fHistList->Add(fh2ESDTriggerVtx);
	211
	212
	213	fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
	214	fHistList->Add(fh1PtHard);
	215	fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
	216	fHistList->Add(fh1PtHardTrials);
f2dd0695	217	fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
f2dd0695	218	AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
bf7b8731	219
f2dd0695	220	fHistList->Add(fh1SelectionInfoESD);
b5a3f310	221	// 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
	222	// 3 triggers BB BE/EB EE
	223
	224	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);
	225	fHistList->Add(fh2ESDTriggerRun);
	226
	227	fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
	228	fHistList->Add(fh2VtxXY);
bf7b8731	229	// =========== Switch on Sumw2 for all histos ===========
	230	for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
	231	TH1 h1 = dynamic_cast<TH1>(fHistList->At(i));
	232	if (h1){
	233	h1->Sumw2();
	234	continue;
	235	}
	236	THnSparse hn = dynamic_cast<THnSparse>(fHistList->At(i));
	237	if(hn)hn->Sumw2();
	238	}
	239
df7848fc	240	fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
	241	fHistList->Add(fh1NCosmicsPerEvent),
	242
bf7b8731	243
	244	TH1::AddDirectory(oldStatus);
	245	}
	246
	247	void AliAnalysisTaskJetServices::Init()
	248	{
	249	//
	250	// Initialization
	251	//
bf7b8731	252	if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
	253
	254	}
	255
	256	void AliAnalysisTaskJetServices::UserExec(Option_t /option*/)
	257	{
	258
	259	//
	260	// Execute analysis for current event
	261	//
	262
	263	AliAODEvent *aod = 0;
	264	AliESDEvent *esd = 0;
	265
	266	AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
f2dd0695	267	fSelectionInfoESD = 0; // reset
	268	AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
	269
bf7b8731	270
	271	if(fUseAODInput){
	272	aod = dynamic_cast<AliAODEvent*>(InputEvent());
	273	if(!aod){
	274	Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
	275	return;
	276	}
	277	// fethc the header
	278	}
	279	else{
	280	// assume that the AOD is in the general output...
	281	aod = AODEvent();
	282	if(!aod){
	283	Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
	284	return;
	285	}
	286	esd = dynamic_cast<AliESDEvent*>(InputEvent());
	287	}
	288
f2dd0695	289	fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kNone;
f2dd0695	290
b5a3f310	291	if(esd){
	292	Float_t run = (Float_t)esd->GetRunNumber();
	293	const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
	294	Float_t zvtx = -999;
	295	Float_t xvtx = -999;
	296	Float_t yvtx = -999;
	297
3568c3d6	298	if(vtxESD->GetNContributors()>2){
b5a3f310	299	zvtx = vtxESD->GetZ();
	300	yvtx = vtxESD->GetY();
	301	xvtx = vtxESD->GetX();
	302	}
	303
	304	Int_t iTrig = -1;
	305	if(esd->GetFiredTriggerClasses().Contains("CINT1B")
	306	\|\|esd->GetFiredTriggerClasses().Contains("CSMBB")
	307	\|\|esd->GetFiredTriggerClasses().Contains("MB1")
	308	\|\|esd->GetFiredTriggerClasses().Contains("CINT6B")){
	309	iTrig = 0;
f2dd0695	310	fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kBunchBunch;
b5a3f310	311	}
	312	else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
	313	\|\|esd->GetFiredTriggerClasses().Contains("CSMBA")
	314	\|\|esd->GetFiredTriggerClasses().Contains("CINT6A")
	315	\|\|esd->GetFiredTriggerClasses().Contains("CINT1C")
	316	\|\|esd->GetFiredTriggerClasses().Contains("CSMBC")
	317	\|\|esd->GetFiredTriggerClasses().Contains("CINT6C")){
	318	// empty bunch or bunch empty
	319	iTrig = 1;
f2dd0695	320	fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kBunchEmpty;
b5a3f310	321	}
f2dd0695	322	else if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
b5a3f310	323	\|\|esd->GetFiredTriggerClasses().Contains("CINT6-E")){
b5a3f310	324	iTrig = 2;
f2dd0695	325	fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kEmptyEmpty;
b5a3f310	326	}
	327
	328
	329	if(iTrig>=0){
	330	iTrig *= 3;
	331	fh2ESDTriggerRun->Fill(run,iTrig+1);
3568c3d6	332	if(vtxESD->GetNContributors()>2){
b5a3f310	333	fh2ESDTriggerRun->Fill(run,iTrig+2);
	334	fh2VtxXY->Fill(xvtx,yvtx);
	335	}
c0997643	336	xvtx -= fVtxXMean;
	337	yvtx -= fVtxYMean;
	338	zvtx -= fVtxZMean;
cf049e94	339	Float_t r2 = xvtx xvtx + yvtx yvtx;
c0997643	340	if(TMath::Abs(zvtx)<fVtxZCut&&r2<fVtxRCut)fh2ESDTriggerRun->Fill(run,iTrig+3);
b5a3f310	341	}
	342	else{
	343	fh2ESDTriggerRun->Fill(run,0);
	344	}
	345
	346
	347	}
bf7b8731	348
bf7b8731	349
c0997643	350	// Apply additional constraints
c0997643	351	Bool_t esdEventSelected = IsEventSelectedESD(esd);
f2dd0695	352	Bool_t esdEventPileUp = IsEventPileUpESD(esd);
	353	Bool_t esdEventCosmic = IsEventCosmicESD(esd);
	354
c0997643	355	Bool_t aodEventSelected = IsEventSelectedAOD(aod);
f2dd0695	356
	357	Bool_t physicsSelection = fInputHandler->IsEventSelected();
	358
	359	if(esdEventSelected) fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kVertexIn;
	360	if(esdEventPileUp) fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kIsPileUp;
	361	if(esdEventCosmic) fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kIsCosmic;
	362	if(physicsSelection) fSelectionInfoESD \|= AliAnalysisHelperJetTasks::kPhysicsSelection;
	363
	364
	365	// here we have all selection information, fill histogram
	366	for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
	367	if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
	368	}
	369	AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
c0997643	370
	371	if(esd&&aod&&fDebug){
	372	if(esdEventSelected&&!aodEventSelected){
	373	Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
	374	const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
bf7b8731	375	const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
c0997643	376	Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
	377	vtxESD->Print();
	378	Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
	379	vtxAOD->Print();
bf7b8731	380	}
c0997643	381	}
	382
	383	// loop over all possible triggers for esd
	384
	385	if(esd){
	386	const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
	387	// Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
	388	TString vtxName(vtxESD->GetName());
	389	for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
	390	Bool_t esdTrig = kFALSE;
	391	esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
	392	if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
216601f0	393	Bool_t cand = fInputHandler->IsEventSelected();
	394	if(cand){
	395	fh2ESDTriggerCount->Fill(it,kSelectedALICE);
216601f0	396	}
cc0649e4	397	if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
3568c3d6	398	if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){
c0997643	399	if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex);
bf7b8731	400	Float_t zvtx = vtxESD->GetZ();
c0997643	401	if(esdEventSelected&&esdTrig){
bf7b8731	402	fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
cc89bb69	403	fh2ESDTriggerVtx->Fill(it,zvtx);
bf7b8731	404	}
bf7b8731	405	}
c0997643	406	if(cand&&esdEventSelected){
	407	fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
	408	fh2ESDTriggerCount->Fill(it,kSelected);
	409	AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
	410	}
bf7b8731	411	}
bf7b8731	412	}
bf7b8731	413
c0997643	414	if(aod){
	415	const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
	416	// Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
	417	TString vtxTitle(vtxAOD->GetTitle());
	418	for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
	419	Bool_t aodTrig = kFALSE;
	420	aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
	421	if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
3568c3d6	422	if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
c0997643	423	if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
c0997643	424	Float_t zvtx = vtxAOD->GetZ();
c0997643	425	if(aodTrig&&aodEventSelected){
cc89bb69	426	fh2TriggerVtx->Fill(it,zvtx);
c0997643	427	fh2TriggerCount->Fill(it,kTriggeredVertexIn);
	428	}
	429	}
	430	}
	431	}
bf7b8731	432
	433	if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
	434
	435
	436	Double_t ptHard = 0;
	437	Double_t nTrials = 1; // Trials for MC trigger
	438
	439	fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
	440	AliMCEvent* mcEvent = MCEvent();
	441	// AliStack *pStack = 0;
	442	if(mcEvent){
	443	AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
	444	if(pythiaGenHeader){
	445	nTrials = pythiaGenHeader->Trials();
	446	ptHard = pythiaGenHeader->GetPtHard();
	447	int iProcessType = pythiaGenHeader->ProcessType();
	448	// 11 f+f -> f+f
	449	// 12 f+barf -> f+barf
	450	// 13 f+barf -> g+g
	451	// 28 f+g -> f+g
	452	// 53 g+g -> f+barf
	453	// 68 g+g -> g+g
	454	if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
	455	if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
	456	fh1PtHard->Fill(ptHard);
	457	fh1PtHardTrials->Fill(ptHard,nTrials);
	458
	459	}// if pythia gen header
	460	}
	461
	462	// trigger selection
	463
	464
	465	PostData(1, fHistList);
	466	}
	467
c0997643	468	Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
	469	if(!esd)return kFALSE;
	470	const AliESDVertex *vtx = esd->GetPrimaryVertex();
	471	// We can check the type of the vertex by its name and title
	472	// if vertexer failed title is empty (default c'tor) and ncontributors is 0
	473	// vtx name title
	474	// Tracks PrimaryVertex VertexerTracksNoConstraint
	475	// TPC TPCVertex VertexerTracksNoConstraint
	476	// SPD SPDVertex vertexer: 3D
	477	// SPD SPDVertex vertexer: Z
	478
	479
	480
3568c3d6	481	if(vtx->GetNContributors()<3)return kFALSE;
c0997643	482
	483	// do not want tpc only primary vertex
	484	TString vtxName(vtx->GetName());
	485	if(vtxName.Contains("TPCVertex"))return kFALSE;
	486	Float_t zvtx = vtx->GetZ();
	487	Float_t yvtx = vtx->GetY();
	488	Float_t xvtx = vtx->GetX();
	489
	490	// here we may fill the histos for run dependence....
	491
	492	xvtx -= fVtxXMean;
	493	yvtx -= fVtxYMean;
	494	zvtx -= fVtxZMean;
	495
	496	Float_t r2 = yvtxyvtx+xvtxxvtx;
	497
	498	Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
	499	return eventSel;
	500	}
	501
f2dd0695	502	Bool_t AliAnalysisTaskJetServices::IsEventPileUpESD(AliESDEvent* esd){
	503	if(!esd)return kFALSE;
	504	return esd->IsPileupFromSPD();
	505	}
	506
	507	Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
	508	if(!esd)return kFALSE;
	509	// add track cuts for which we look for cosmics...
df7848fc	510
	511	Bool_t isCosmic = kFALSE;
	512	Int_t nTracks = esd->GetNumberOfTracks();
	513	Int_t nCosmicCandidates = 0;
	514
	515	for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
	516	AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
	517	if (!track1) continue;
	518	if(track1->Pt()<fPtMinCosmic) continue;
	519	//Start 2nd track loop to look for correlations
	520	for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
	521	AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
	522	if(!track2) continue;
	523	if(track2->Pt()<fPtMinCosmic) continue;
	524	//Check if back-to-back
	525	Double_t mom1[3],mom2[3];
	526	track1->GetPxPyPz(mom1);
	527	track2->GetPxPyPz(mom2);
	528	TVector3 momv1(mom1[0],mom1[1],mom1[2]);
	529	TVector3 momv2(mom2[0],mom2[1],mom2[2]);
	530	Float_t theta = (float)(momv1.Phi()-momv2.Phi());
	531	if(theta<-0.5TMath::Pi()) theta+=2.TMath::Pi();
	532
	533	Float_t deltaPhi = track1->Phi()-track2->Phi();
	534	if(deltaPhi<-0.5TMath::Pi()) deltaPhi+=2.TMath::Pi();
	535
	536	Float_t rIsol = (float)(TMath::Sqrt( deltaPhideltaPhi+(track1->Eta()-track2->Eta())(track1->Eta()-track2->Eta()) ));
	537	if(rIsol<fRIsolMinCosmic) continue;
	538
	539	if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
	540	nCosmicCandidates+=1.;
	541	isCosmic = kTRUE;
	542	}
	543
	544	}
	545	}
	546
	547	fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
	548
	549	return isCosmic;
f2dd0695	550	}
	551
	552
c0997643	553	Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
	554	if(!aod)return kFALSE;
	555	const AliAODVertex *vtx = aod->GetPrimaryVertex();
	556	// We can check the type of the vertex by its name and title
	557	// if vertexer failed title is empty (default c'tor) and ncontributors is 0
	558	// vtx name title
	559	// Tracks PrimaryVertex VertexerTracksNoConstraint
	560	// TPC TPCVertex VertexerTracksNoConstraint
	561	// SPD SPDVertex vertexer: 3D
	562	// SPD SPDVertex vertexer: Z
	563
	564
	565
3568c3d6	566	if(vtx->GetNContributors()<3)return kFALSE;
c0997643	567
	568	// do not want tpc only primary vertex
	569	TString vtxName(vtx->GetName());
	570	if(vtxName.Contains("TPCVertex"))return kFALSE;
	571
	572	Float_t zvtx = vtx->GetZ();
	573	Float_t yvtx = vtx->GetY();
	574	Float_t xvtx = vtx->GetX();
	575
	576	// here we may fill the histos for run dependence....
	577
	578	xvtx -= fVtxXMean;
	579	yvtx -= fVtxYMean;
	580	zvtx -= fVtxZMean;
	581
	582	Float_t r2 = yvtxyvtx+xvtxxvtx;
	583	Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
	584	return eventSel;
	585	}
	586
	587
bf7b8731	588
f2dd0695	589
bf7b8731	590	void AliAnalysisTaskJetServices::Terminate(Option_t /option*/)
	591	{
	592	// Terminate analysis
	593	//
bf7b8731	594	}