[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),
  fAvgTrials(1),
  fVtxXMean(0),
  fVtxYMean(0),
  fVtxZMean(0),
  fVtxRCut(1.),
  fVtxZCut(8.),
  fh1Xsec(0x0),
  fh1Trials(0x0),
  fh1PtHard(0x0),
  fh1PtHardTrials(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),
  fAvgTrials(1),
  fVtxXMean(0),
  fVtxYMean(0),
  fVtxZMean(0),
  fVtxRCut(1.),
  fVtxZCut(8.),
  fh1Xsec(0x0),
  fh1Trials(0x0),
  fh1PtHard(0x0),
  fh1PtHardTrials(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);

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

  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());
  }
  
  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;
    }
    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;
    }
    if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
       ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
      iTrig = 2;
    }

    
    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 aodEventSelected = IsEventSelectedAOD(aod);
  

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