[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 <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),
  fZVtxCut(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),
  fZVtxCut(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
  //

  Printf(">>> AnalysisTaskJetServices::Init() debug level %d\n",fDebug);
  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()>0){
      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()>0){
	fh2ESDTriggerRun->Fill(run,iTrig+2);
	fh2VtxXY->Fill(xvtx,yvtx);
      }
      Float_t r2 = xvtx *xvtx + yvtx *yvtx; 
      if(TMath::Abs(zvtx)<fZVtxCut&&r2<1)fh2ESDTriggerRun->Fill(run,iTrig+3);
    }
    else{
      fh2ESDTriggerRun->Fill(run,0);
    }


  }
  
  // loop over all possible trigger and 
  for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
    Bool_t esdTrig = kFALSE;
    if(esd){
      esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
      if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
    }
    Bool_t aodTrig = kFALSE;
    if(aod){
      aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
      if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
    }

    // Check wether we have also an SPD vertex
    
    if(aod){
      const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
      //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
      
      if(vtxAOD->GetNContributors()>0){
	if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
	Float_t zvtx = vtxAOD->GetZ();
	Float_t yvtx = vtxAOD->GetY();
	Float_t xvtx = vtxAOD->GetX();
	fh2TriggerVtx->Fill(it,zvtx);
	Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
	if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&r2<1){
	  fh2TriggerCount->Fill(it,kTriggeredVertexIn);
	}
      }
    }
    if(esd){
      const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
      //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
      Bool_t cand = fInputHandler->IsEventSelected();
      if(cand){
	fh2ESDTriggerCount->Fill(it,kSelectedALICE); 
      }
      if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);

      if(vtxESD->GetNContributors()>0){
	if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
	Float_t zvtx = vtxESD->GetZ();
    	Float_t yvtx = vtxESD->GetY();
    	Float_t xvtx = vtxESD->GetX();
	Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
	fh2ESDTriggerVtx->Fill(it,zvtx);
	if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&r2<1){
	  fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
	  if(cand){
	    fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
	    fh2ESDTriggerCount->Fill(it,kSelected);
	    AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
	  }
	  // here we select based on ESD info...
	}
      }
    }
  }


  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);
}


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>
	25	#include <TSystem.h>
	26	#include <TInterpreter.h>
	27	#include <TChain.h>
	28	#include <TFile.h>
	29	#include <TKey.h>
	30	#include <TH1F.h>
	31	#include <TH2F.h>
	32	#include <TH3F.h>
	33	#include <TProfile.h>
	34	#include <TList.h>
	35	#include <TLorentzVector.h>
	36	#include <TClonesArray.h>
	37	#include "TDatabasePDG.h"
	38
	39	#include "AliAnalysisTaskJetServices.h"
fe669ac6	40	#include "AliAnalysisDataContainer.h"
fe669ac6	41	#include "AliAnalysisDataSlot.h"
bf7b8731	42	#include "AliAnalysisManager.h"
	43	#include "AliJetFinder.h"
	44	#include "AliJetHeader.h"
	45	#include "AliJetReader.h"
	46	#include "AliJetReaderHeader.h"
	47	#include "AliUA1JetHeaderV1.h"
	48	#include "AliESDEvent.h"
	49	#include "AliAODEvent.h"
	50	#include "AliAODHandler.h"
	51	#include "AliAODTrack.h"
	52	#include "AliAODJet.h"
	53	#include "AliAODMCParticle.h"
	54	#include "AliMCEventHandler.h"
	55	#include "AliMCEvent.h"
	56	#include "AliStack.h"
	57	#include "AliGenPythiaEventHeader.h"
	58	#include "AliJetKineReaderHeader.h"
	59	#include "AliGenCocktailEventHeader.h"
	60	#include "AliInputEventHandler.h"
fe669ac6	61	#include "AliPhysicsSelection.h"
bf7b8731	62
	63
	64	#include "AliAnalysisHelperJetTasks.h"
	65
	66	ClassImp(AliAnalysisTaskJetServices)
	67
	68	AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
	69	fUseAODInput(kFALSE),
cc0649e4	70	fUsePhysicsSelection(kFALSE),
cc0649e4	71	fRealData(kFALSE),
bf7b8731	72	fAvgTrials(1),
fe669ac6	73	fZVtxCut(8.),
bf7b8731	74	fh1Xsec(0x0),
	75	fh1Trials(0x0),
	76	fh1PtHard(0x0),
	77	fh1PtHardTrials(0x0),
	78	fh2TriggerCount(0x0),
	79	fh2ESDTriggerCount(0x0),
	80	fh2TriggerVtx(0x0),
	81	fh2ESDTriggerVtx(0x0),
b5a3f310	82	fh2ESDTriggerRun(0x0),
b5a3f310	83	fh2VtxXY(0x0),
bf7b8731	84	fHistList(0x0)
bf7b8731	85	{
b5a3f310	86	fRunRange[0] = fRunRange[1] = 0;
bf7b8731	87	}
	88
	89	AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
	90	AliAnalysisTaskSE(name),
	91	fUseAODInput(kFALSE),
cc0649e4	92	fUsePhysicsSelection(kFALSE),
cc0649e4	93	fRealData(kFALSE),
bf7b8731	94	fAvgTrials(1),
fe669ac6	95	fZVtxCut(8.),
bf7b8731	96	fh1Xsec(0x0),
	97	fh1Trials(0x0),
	98	fh1PtHard(0x0),
	99	fh1PtHardTrials(0x0),
	100	fh2TriggerCount(0x0),
	101	fh2ESDTriggerCount(0x0),
	102	fh2TriggerVtx(0x0),
	103	fh2ESDTriggerVtx(0x0),
b5a3f310	104	fh2ESDTriggerRun(0x0),
b5a3f310	105	fh2VtxXY(0x0),
bf7b8731	106	fHistList(0x0)
bf7b8731	107	{
b5a3f310	108	fRunRange[0] = fRunRange[1] = 0;
bf7b8731	109	DefineOutput(1,TList::Class());
	110	}
	111
	112
	113
	114	Bool_t AliAnalysisTaskJetServices::Notify()
	115	{
	116	//
	117	// Implemented Notify() to read the cross sections
	118	// and number of trials from pyxsec.root
	119	//
	120
	121	TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
	122	Float_t xsection = 0;
	123	Float_t ftrials = 1;
	124
	125	fAvgTrials = 1;
	126	if(tree){
	127	TFile *curfile = tree->GetCurrentFile();
	128	if (!curfile) {
	129	Error("Notify","No current file");
	130	return kFALSE;
	131	}
	132	if(!fh1Xsec\|\|!fh1Trials){
	133	Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
	134	return kFALSE;
	135	}
	136	AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
	137	fh1Xsec->Fill("<#sigma>",xsection);
	138	// construct a poor man average trials
	139	Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
	140	if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
	141	}
	142	return kTRUE;
	143	}
	144
	145	void AliAnalysisTaskJetServices::UserCreateOutputObjects()
	146	{
	147
	148	//
	149	// Create the output container
	150	//
	151
	152
bf7b8731	153	if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
	154
	155	OpenFile(1);
	156	if(!fHistList)fHistList = new TList();
	157
	158	Bool_t oldStatus = TH1::AddDirectoryStatus();
	159	TH1::AddDirectory(kFALSE);
bf7b8731	160	fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
	161	fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
	162	fHistList->Add(fh1Xsec);
	163
	164	fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
	165	fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
	166	fHistList->Add(fh1Trials);
	167
	168	const Int_t nBinPt = 100;
	169	Double_t binLimitsPt[nBinPt+1];
	170	for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
	171	if(iPt == 0){
	172	binLimitsPt[iPt] = 0.0;
	173	}
	174	else {// 1.0
	175	binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
	176	}
	177	}
	178
	179	fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
	180	fHistList->Add(fh2TriggerCount);
	181
	182	fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
	183	fHistList->Add(fh2ESDTriggerCount);
	184
	185	fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
	186	fHistList->Add(fh2TriggerVtx);
	187
	188	fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
	189	fHistList->Add(fh2ESDTriggerVtx);
	190
	191
	192	fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
	193	fHistList->Add(fh1PtHard);
	194	fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
	195	fHistList->Add(fh1PtHardTrials);
	196
b5a3f310	197	// 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
	198	// 3 triggers BB BE/EB EE
	199
	200	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);
	201	fHistList->Add(fh2ESDTriggerRun);
	202
	203	fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
	204	fHistList->Add(fh2VtxXY);
bf7b8731	205	// =========== Switch on Sumw2 for all histos ===========
	206	for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
	207	TH1 h1 = dynamic_cast<TH1>(fHistList->At(i));
	208	if (h1){
	209	h1->Sumw2();
	210	continue;
	211	}
	212	THnSparse hn = dynamic_cast<THnSparse>(fHistList->At(i));
	213	if(hn)hn->Sumw2();
	214	}
	215
	216
	217	TH1::AddDirectory(oldStatus);
	218	}
	219
	220	void AliAnalysisTaskJetServices::Init()
	221	{
	222	//
	223	// Initialization
	224	//
	225
	226	Printf(">>> AnalysisTaskJetServices::Init() debug level %d\n",fDebug);
	227	if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
	228
	229	}
	230
	231	void AliAnalysisTaskJetServices::UserExec(Option_t /option*/)
	232	{
	233
	234	//
	235	// Execute analysis for current event
	236	//
	237
	238	AliAODEvent *aod = 0;
	239	AliESDEvent *esd = 0;
	240
	241	AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
	242
	243	if(fUseAODInput){
	244	aod = dynamic_cast<AliAODEvent*>(InputEvent());
	245	if(!aod){
	246	Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
	247	return;
	248	}
	249	// fethc the header
	250	}
	251	else{
	252	// assume that the AOD is in the general output...
	253	aod = AODEvent();
	254	if(!aod){
	255	Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
	256	return;
	257	}
	258	esd = dynamic_cast<AliESDEvent*>(InputEvent());
	259	}
	260
b5a3f310	261	if(esd){
	262	Float_t run = (Float_t)esd->GetRunNumber();
	263	const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
	264	Float_t zvtx = -999;
	265	Float_t xvtx = -999;
	266	Float_t yvtx = -999;
	267
	268	if(vtxESD->GetNContributors()>0){
	269	zvtx = vtxESD->GetZ();
	270	yvtx = vtxESD->GetY();
	271	xvtx = vtxESD->GetX();
	272	}
	273
	274	Int_t iTrig = -1;
	275	if(esd->GetFiredTriggerClasses().Contains("CINT1B")
	276	\|\|esd->GetFiredTriggerClasses().Contains("CSMBB")
	277	\|\|esd->GetFiredTriggerClasses().Contains("MB1")
	278	\|\|esd->GetFiredTriggerClasses().Contains("CINT6B")){
	279	iTrig = 0;
	280	}
	281	else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
	282	\|\|esd->GetFiredTriggerClasses().Contains("CSMBA")
	283	\|\|esd->GetFiredTriggerClasses().Contains("CINT6A")
	284	\|\|esd->GetFiredTriggerClasses().Contains("CINT1C")
	285	\|\|esd->GetFiredTriggerClasses().Contains("CSMBC")
	286	\|\|esd->GetFiredTriggerClasses().Contains("CINT6C")){
	287	// empty bunch or bunch empty
	288	iTrig = 1;
	289	}
	290	if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
	291	\|\|esd->GetFiredTriggerClasses().Contains("CINT6-E")){
	292	iTrig = 2;
	293	}
	294
	295
	296	if(iTrig>=0){
	297	iTrig *= 3;
	298	fh2ESDTriggerRun->Fill(run,iTrig+1);
	299	if(vtxESD->GetNContributors()>0){
	300	fh2ESDTriggerRun->Fill(run,iTrig+2);
	301	fh2VtxXY->Fill(xvtx,yvtx);
	302	}
cf049e94	303	Float_t r2 = xvtx xvtx + yvtx yvtx;
cf049e94	304	if(TMath::Abs(zvtx)<fZVtxCut&&r2<1)fh2ESDTriggerRun->Fill(run,iTrig+3);
b5a3f310	305	}
	306	else{
	307	fh2ESDTriggerRun->Fill(run,0);
	308	}
	309
	310
	311	}
bf7b8731	312
	313	// loop over all possible trigger and
	314	for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
	315	Bool_t esdTrig = kFALSE;
	316	if(esd){
	317	esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
	318	if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
	319	}
	320	Bool_t aodTrig = kFALSE;
	321	if(aod){
	322	aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
	323	if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
	324	}
	325
	326	// Check wether we have also an SPD vertex
	327
	328	if(aod){
	329	const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
	330	// Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
	331
	332	if(vtxAOD->GetNContributors()>0){
	333	if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
	334	Float_t zvtx = vtxAOD->GetZ();
	335	Float_t yvtx = vtxAOD->GetY();
	336	Float_t xvtx = vtxAOD->GetX();
	337	fh2TriggerVtx->Fill(it,zvtx);
cf049e94	338	Float_t r2 = yvtxyvtx+xvtxxvtx;
cf049e94	339	if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&r2<1){
bf7b8731	340	fh2TriggerCount->Fill(it,kTriggeredVertexIn);
	341	}
	342	}
	343	}
	344	if(esd){
	345	const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
	346	// Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
216601f0	347	Bool_t cand = fInputHandler->IsEventSelected();
	348	if(cand){
	349	fh2ESDTriggerCount->Fill(it,kSelectedALICE);
216601f0	350	}
cc0649e4	351	if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
cc0649e4	352
bf7b8731	353	if(vtxESD->GetNContributors()>0){
	354	if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
	355	Float_t zvtx = vtxESD->GetZ();
	356	Float_t yvtx = vtxESD->GetY();
	357	Float_t xvtx = vtxESD->GetX();
a923bd34	358	Float_t r2 = yvtxyvtx+xvtxxvtx;
bf7b8731	359	fh2ESDTriggerVtx->Fill(it,zvtx);
a923bd34	360	if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&r2<1){
bf7b8731	361	fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
a923bd34	362	if(cand){
	363	fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
	364	fh2ESDTriggerCount->Fill(it,kSelected);
	365	AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
	366	}
bf7b8731	367	// here we select based on ESD info...
bf7b8731	368	}
bf7b8731	369	}
bf7b8731	370	}
bf7b8731	371	}
	372
	373
	374
	375	if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
	376
	377
	378	Double_t ptHard = 0;
	379	Double_t nTrials = 1; // Trials for MC trigger
	380
	381	fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
	382	AliMCEvent* mcEvent = MCEvent();
	383	// AliStack *pStack = 0;
	384	if(mcEvent){
	385	AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
	386	if(pythiaGenHeader){
	387	nTrials = pythiaGenHeader->Trials();
	388	ptHard = pythiaGenHeader->GetPtHard();
	389	int iProcessType = pythiaGenHeader->ProcessType();
	390	// 11 f+f -> f+f
	391	// 12 f+barf -> f+barf
	392	// 13 f+barf -> g+g
	393	// 28 f+g -> f+g
	394	// 53 g+g -> f+barf
	395	// 68 g+g -> g+g
	396	if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
	397	if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
	398	fh1PtHard->Fill(ptHard);
	399	fh1PtHardTrials->Fill(ptHard,nTrials);
	400
	401	}// if pythia gen header
	402	}
	403
	404	// trigger selection
	405
	406
	407	PostData(1, fHistList);
	408	}
	409
	410
	411	void AliAnalysisTaskJetServices::Terminate(Option_t /option*/)
	412	{
	413	// Terminate analysis
	414	//
bf7b8731	415	}