+++ /dev/null
-#ifndef ALIANALYSISTASKSE_H
-
-#include <Riostream.h>
-#include <TROOT.h>
-#include <TFile.h>
-#include <TChain.h>
-#include <TTree.h>
-#include <TKey.h>
-#include <TProfile.h>
-#include <TProfile2D.h>
-#include <TH1.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TH1D.h>
-#include <TH1I.h>
-#include <TArrayF.h>
-#include <THnSparse.h>
-#include <TCanvas.h>
-#include <TList.h>
-#include <TClonesArray.h>
-#include <TObject.h>
-#include <TMath.h>
-#include <TSystem.h>
-#include <TInterpreter.h>
-#include "AliAnalysisTask.h"
-#include "AliCentrality.h"
-#include "AliStack.h"
-#include "AliESDEvent.h"
-#include "AliESDInputHandler.h"
-#include "AliAODEvent.h"
-#include "AliAODHandler.h"
-#include "AliAnalysisManager.h"
-#include "AliAnalysisTaskSE.h"
-#include "AliAnalysisHelperJetTasks.h"
-#include "AliInputEventHandler.h"
-#endif
-
-#include <time.h>
-#include <TRandom3.h>
-#include "AliGenEventHeader.h"
-#include "AliGenPythiaEventHeader.h"
-#include "AliGenHijingEventHeader.h"
-#include "AliAODMCHeader.h"
-#include "AliMCEvent.h"
-#include "AliLog.h"
-#include <AliEmcalJet.h>
-#include <AliPicoTrack.h>
-#include "AliVEventHandler.h"
-#include "AliVParticle.h"
-#include "AliAODMCParticle.h"
-#include "AliAnalysisUtils.h"
-#include "AliRhoParameter.h"
-#include "TVector3.h"
-#include "AliVVertex.h"
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <vector>
-
-#include "AliAnalysisTaskHJetSpectra.h"
-using std::min;
-
-// ANALYSIS OF HIGH PT HADRON TRIGGER ASSOCIATED SPECTRUM OF RECOIL JETS IN P+PB
-// Author Filip Krizek (17.May. 2014)
-
-//TODO: Not accessing the particles when using MC
-//TODO: FillHistogram can be done better with virtual TH1(?)
-ClassImp(AliAnalysisTaskHJetSpectra)
-//________________________________________________________________________________________
-
-AliAnalysisTaskHJetSpectra::AliAnalysisTaskHJetSpectra():
-AliAnalysisTaskSE(), fOutputList(0), fAnalyzePythia(0), fAnalyzeHijing(0), fIsKinematics(0), fUseDefaultVertexCut(1), fUsePileUpCut(1),
-fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fRhoTaskName(),
-fRandConeRadius(0.4),fRandConeRadiusSquared(fRandConeRadius*fRandConeRadius), fSignalJetRadius(0.4), fBackgroundJetRadius(0.3), fBackgroundJetPtMin(15.0),
-fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.5), fNumberOfCentralityBins(20), fCentralityType("V0A"),
-fCrossSection(0.0), fTrials(0.0), fImpParam(-1.0), fRandom(0), fHelperClass(0), fInitialized(0),
-fTTlow(8.0), fTThigh(9.0), fTTtype(0), fDphiCut(TMath::Pi()-0.6), fUseDoubleBinPrecision(0),
-fHistEvtSelection(0x0), fh2Ntriggers(0x0), fHJetSpec(0x0), fHJetSpecSubUeMedian(0x0), fHJetSpecSubUeCone(0x0), fHJetSpecSubUeCMS(0x0),
-fhRhoCellMedian(0x0), fhRhoCone(0x0), fhRhoCMS(0x0),
-fhRhoCellMedianIncl(0x0), fhRhoConeIncl(0x0), fhRhoCMSIncl(0x0),
-fARhoCellMedian(0x0), fARhoCone(0x0), fARhoCMS(0x0),
-fhDeltaPtMedian(0x0), fhDeltaPtCone(0x0), fhDeltaPtCMS(0x0),
-fhDeltaPtMedianIncl(0x0), fhDeltaPtConeIncl(0x0), fhDeltaPtCMSIncl(0x0),
-fhDeltaPtMedianNearSide(0x0), fhDeltaPtMedianAwaySide(0x0), fhDeltaPtCMSNearSide(0x0), fhDeltaPtCMSAwaySide(0x0),
-fhDeltaPtMedianExclTrigCone(0x0),fhDeltaPtCMSExclTrigCone(0x0), fhDeltaPtMedianExclAwayJet(0x0), fhDeltaPtCMSExclAwayJet(0x0),
-fhJetPhi(0x0), fhTrackPhi(0x0), fhJetEta(0x0), fhTrackEta(0x0), fhTrackCentVsPt(0x0), fhVertexZ(0x0), fhVertexZAccept(0x0),
-fhDphiTriggerJetMinBias(0x0),fhDphiTriggerJetCent20(0x0), fhDphiTriggerJetAccept(0x0),
-fhCentrality(0x0), fhCentralityV0M(0x0), fhCentralityV0A(0x0), fhCentralityV0C(0x0), fhCentralityZNA(0x0),
-fNofRndTrials(2000), fJetFreeAreaFrac(0.8), fnEta(2), fnPhi(11), fEtaSize(0.9), fPhiSize(2*TMath::Pi()/fnPhi), fCellArea(fPhiSize*fEtaSize),
-fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), fhImpactParameter(0x0), fhImpactParameterTT(0x0),
-fNofRandomCones(1),
-fRConesR(0.1),fRConesRSquared(fRConesR*fRConesR),fnRCones(16)
-{
- //default constructor
- for(Int_t k=0; k<50; k++){
- fRConePhi[k] = 0.0;
- fRConeEta[k] = 0.0;
- }
-}
-
-//________________________________________________________________________
-AliAnalysisTaskHJetSpectra::AliAnalysisTaskHJetSpectra(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) :
-AliAnalysisTaskSE(name), fOutputList(0), fAnalyzePythia(0), fAnalyzeHijing(0), fIsKinematics(0), fUseDefaultVertexCut(1), fUsePileUpCut(1),
- fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fRhoTaskName(),
-fRandConeRadius(0.4),fRandConeRadiusSquared(fRandConeRadius*fRandConeRadius), fSignalJetRadius(0.4), fBackgroundJetRadius(0.3), fBackgroundJetPtMin(15.0),
-fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.5), fNumberOfCentralityBins(20), fCentralityType("V0A"),
-fCrossSection(0.0), fTrials(0.0), fImpParam(-1.0), fRandom(0), fHelperClass(0), fInitialized(0),
-fTTlow(8.0), fTThigh(9.0), fTTtype(0), fDphiCut(TMath::Pi()-0.6), fUseDoubleBinPrecision(0),
-fHistEvtSelection(0x0), fh2Ntriggers(0x0), fHJetSpec(0x0), fHJetSpecSubUeMedian(0x0), fHJetSpecSubUeCone(0x0), fHJetSpecSubUeCMS(0x0),
-fhRhoCellMedian(0x0), fhRhoCone(0x0), fhRhoCMS(0x0),
-fhRhoCellMedianIncl(0x0), fhRhoConeIncl(0x0), fhRhoCMSIncl(0x0),
-fARhoCellMedian(0x0), fARhoCone(0x0), fARhoCMS(0x0),
-fhDeltaPtMedian(0x0), fhDeltaPtCone(0x0), fhDeltaPtCMS(0x0),
-fhDeltaPtMedianIncl(0x0), fhDeltaPtConeIncl(0x0), fhDeltaPtCMSIncl(0x0),
-fhDeltaPtMedianNearSide(0x0), fhDeltaPtMedianAwaySide(0x0), fhDeltaPtCMSNearSide(0x0), fhDeltaPtCMSAwaySide(0x0),
-fhDeltaPtMedianExclTrigCone(0x0),fhDeltaPtCMSExclTrigCone(0x0), fhDeltaPtMedianExclAwayJet(0x0), fhDeltaPtCMSExclAwayJet(0x0),
-fhJetPhi(0x0), fhTrackPhi(0x0), fhJetEta(0x0), fhTrackEta(0x0), fhTrackCentVsPt(0x0), fhVertexZ(0x0), fhVertexZAccept(0x0),
-fhDphiTriggerJetMinBias(0x0), fhDphiTriggerJetCent20(0x0), fhDphiTriggerJetAccept(0x0),
-fhCentrality(0x0), fhCentralityV0M(0x0), fhCentralityV0A(0x0), fhCentralityV0C(0x0), fhCentralityZNA(0x0),
-fNofRndTrials(2000), fJetFreeAreaFrac(0.8), fnEta(2), fnPhi(11), fEtaSize(0.9), fPhiSize(2*TMath::Pi()/fnPhi), fCellArea(fPhiSize*fEtaSize),
-fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), fhImpactParameter(0x0), fhImpactParameterTT(0x0),
-fNofRandomCones(1),
-fRConesR(0.1),fRConesRSquared(fRConesR*fRConesR), fnRCones(16)
-{
- //constructor that is called
- //LIST OF TRACKS
- fTrackArrayName = new TString(trackArrayName);
- if((fTrackArrayName->Contains("MC") && fTrackArrayName->Contains("Particles")) ||
- (fTrackArrayName->Contains("mc") && fTrackArrayName->Contains("particles"))){
- fIsKinematics = kTRUE;
- }
-
- //LIST of JETS
- fJetArrayName = new TString(jetArrayName);
- if(strcmp(fJetArrayName->Data(),"") == 0){
- AliError(Form("%s: Jet branch missing !", GetName()));
- }
-
- //LIST OF JETS TO BE IGNORED WHILE RHO ESTIMATE
- fBackgroundJetArrayName = new TString(backgroundJetArrayName); //jets to be removed from cell median rho estimate
- if(strcmp(fBackgroundJetArrayName->Data(),"") == 0){
- AliError(Form("%s: Bg Jet branch missing !", GetName()));
- }
-
- for(Int_t k=0; k<50; k++){
- fRConePhi[k] = 0.0;
- fRConeEta[k] = 0.0;
- }
-
- DefineOutput(1, TList::Class());
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskHJetSpectra::SetAnalyzeMC(Int_t val){
- if(val==1){
- fAnalyzePythia = kTRUE;
- fAnalyzeHijing = kFALSE;
- return;
- }
- if(val==2){
- fAnalyzeHijing = kTRUE;
- fAnalyzePythia = kFALSE;
- return;
- }
-
- fAnalyzeHijing = kFALSE;
- fAnalyzePythia = kFALSE;
- return;
-}
-//________________________________________________________________________
-Double_t AliAnalysisTaskHJetSpectra::GetConePt(Double_t eta, Double_t phi, Double_t radius){
- //sum up pt inside a cone
- Double_t tmpConePt = 0.0;
- Double_t dphi = 0.0;
- Double_t deta = 0.0;
- Double_t radiussquared = radius*radius;
-
- for(Int_t i = 0; i < fTrackArray->GetEntries(); i++){
- AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
- if(!tmpTrack) continue;
- if(IsTrackInAcceptance(tmpTrack)){
- dphi = RelativePhi(tmpTrack->Phi(),phi);
- deta = tmpTrack->Eta() - eta;
- if( dphi*dphi + deta*deta < radiussquared ){
- tmpConePt = tmpConePt + tmpTrack->Pt();
- }
- }
- }
- return tmpConePt;
-}
-
-
-//________________________________________________________________________
-Double_t AliAnalysisTaskHJetSpectra::GetPtHard(){
- //get pt hard from pythia
- AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
- if(MCEvent()){
- if(!pythiaHeader){
- // Check if AOD
- AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
-
- if(aodMCH){
- for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++){
- pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
- if(pythiaHeader) break;
- }
- }
- }
- }
- if(pythiaHeader){
-
- TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
-
- if(tree){
- TFile *curfile = tree->GetCurrentFile();
- if(!curfile) {
- Error("Notify","No current file");
- return kFALSE;
- }
- Float_t xsection = 0.0;
- Float_t trials = 1.0;
-
- AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
-
- fCrossSection = (Double_t) xsection;//pythiaHeader->GetXsection();
-
- if(fCrossSection>0.){ //save cross-section and the number of trials
- fTrials = (Double_t) trials; //pythiaHeader->Trials();
- fh1Xsec->Fill("<#sigma>", fCrossSection);
- fh1Trials->Fill("#sum{ntrials}",fTrials);
- }
- }
- return pythiaHeader->GetPtHard();
- }
- AliWarning(Form("In task %s: GetPtHard() failed!", GetName()));
- return -1.0;
-}
-//________________________________________________________________________
-Double_t AliAnalysisTaskHJetSpectra::GetImpactParameter(){
- //get impact parameter from hijing
- AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader());
- if(MCEvent()){
- if(!hijingHeader){
- // Check if AOD
- AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
-
- if(aodMCH){
- for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++){
- hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(aodMCH->GetCocktailHeader(i));
- if(hijingHeader) break;
- }
- }
- }
- }
- if(hijingHeader){
- return (Double_t) (hijingHeader->ImpactParameter());
- }
- AliWarning(Form("In task %s: GetImpactParameter() failed!", GetName()));
- return -1.0;
-}
-//________________________________________________________________________
-Double_t AliAnalysisTaskHJetSpectra::GetSimPrimaryVertex(){
- //get generator level primary vertex
- AliGenEventHeader* mcHeader = NULL;
- AliAODMCHeader* aodMCH = NULL;
- if(MCEvent()){
- if(fAnalyzePythia){
- mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
- if(!mcHeader){
- // Check if AOD
- aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
-
- if(aodMCH){
- for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){
- mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
- if(mcHeader) break;
- }
- }
- }
- }
-
- if(fAnalyzeHijing){
- mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader());
- if(!mcHeader){
- // Check if AOD
- aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
-
- if(aodMCH){
- for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){
- mcHeader = dynamic_cast<AliGenHijingEventHeader*>(aodMCH->GetCocktailHeader(i));
- if(mcHeader) break;
- }
- }
- }
- }
- }
- if(mcHeader){
-
- TArrayF pyVtx(3);
- mcHeader->PrimaryVertex(pyVtx);
- return (Double_t) (pyVtx[2]);
- }
- AliWarning(Form("In task %s: Pythia Vertex failed!", GetName()));
- return 9999.0;
-}
-
-
-
-//________________________________________________________________________
-/*Double_t AliAnalysisTaskHJetSpectra::GetPythiaTrials()
-{
- #ifdef DEBUGMODE
- AliInfo("Starting GetPythiaTrials.");
- #endif
- AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
- if (MCEvent())
- if (!pythiaHeader)
- {
- // Check if AOD
- AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
-
- if (aodMCH)
- {
- for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++)
- {
- pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
- if (pythiaHeader) break;
- }
- }
- }
-
- #ifdef DEBUGMODE
- AliInfo("Ending GetPythiaTrials.");
- #endif
- if (pythiaHeader)
- return pythiaHeader->Trials();
-
- AliWarning(Form("In task %s: GetPythiaTrials() failed!", GetName()));
- return -1.0;
-}
-*/
-
-//________________________________________________________________________
-Double_t AliAnalysisTaskHJetSpectra::GetExternalRho(){
- // Get rho from event using CMS approach
-
- AliRhoParameter *rho = 0;
- if(!fRhoTaskName.IsNull()) {
- rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoTaskName.Data()));
- if(!rho){
- AliWarning(Form("%s: Could not retrieve rho with name %s!", GetName(), fRhoTaskName.Data()));
- return 0.0;
- }
- }else return 0.0;
-
- return (rho->GetVal());
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskHJetSpectra::IsEventInAcceptance(AliVEvent* event){
- //EVENT SELECTION
-
-
- if(!event) return kFALSE;
-
- //___________________________________________________
-
- if(fAnalyzePythia || fAnalyzeHijing){ //PURE MC
- if(!MCEvent()) return kFALSE;
-
- //BEFORE VERTEX CUT
- Double_t vtxMC = GetSimPrimaryVertex();
- fhVertexZ->Fill(vtxMC);
-
- if(TMath::Abs(vtxMC) > 10.0){
- fHistEvtSelection->Fill(3); //count events rejected by vertex cut
- return kFALSE;
- }
- fhVertexZAccept->Fill(vtxMC);
-
- return kTRUE;
- }
- //___________________________________________________
- //TEST PILE UP
- if(fUsePileUpCut){
- if(!fHelperClass || fHelperClass->IsPileUpEvent(event)){
- fHistEvtSelection->Fill(2); //count events rejected by pileup
- return kFALSE;
- }
- }
- //___________________________________________________
- //BEFORE VERTEX CUT
- fhVertexZ->Fill(event->GetPrimaryVertex()->GetZ());
-
- if(fUseDefaultVertexCut){
- if(!fHelperClass || !fHelperClass->IsVertexSelected2013pA(event)){
- fHistEvtSelection->Fill(3); //count events rejected by vertex cut
- return kFALSE;
- }
- }else{
- if(TMath::Abs(event->GetPrimaryVertex()->GetZ()) > 10.0){
- fHistEvtSelection->Fill(3); //count events rejected by vertex cut
- return kFALSE;
- }
- }
- //___________________________________________________
- //AFTER VERTEX CUT
- fhVertexZAccept->Fill(event->GetPrimaryVertex()->GetZ());
-
- return kTRUE;
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskHJetSpectra::IsTrackInAcceptance(AliVParticle* track){
- // Check if the track pt and eta range
- if(track != 0){
- if(fIsKinematics){
- // TODO: Only working for AOD MC
- if((!track->Charge()) || (!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) )
- return kFALSE;
- }
- if(TMath::Abs(track->Eta()) <= fTrackEtaWindow){ //APPLY TRACK ETA CUT
- if(track->Pt() >= fMinTrackPt){ //APPLY TRACK CUT
- return kTRUE;
- }
- }
- }
- return kFALSE;
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskHJetSpectra::IsBackgroundJetInAcceptance(AliEmcalJet *jet){
- //find jets to be removed from bg calculation
- if(jet != 0){
- if(TMath::Abs(jet->Eta()) <= fBackgroundJetEtaWindow){
- if(jet->Pt() >= fBackgroundJetPtMin){ //accept only hard jets
- return kTRUE;
- }
- }
- }
- return kFALSE;
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskHJetSpectra::IsSignalJetInAcceptance(AliEmcalJet *jet){
- //select jets in acceptance
- if(jet == 0) return kFALSE;
- if(TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow){
- if(jet->Pt() >= fMinTrackPt){
- if(jet->Area() >= fMinJetArea){
- return kTRUE;
- }
- }
- }
- return kFALSE;
-}
-
-
-//________________________________________________________________________
-void AliAnalysisTaskHJetSpectra::ExecOnce(){
- //Read arrays of jets and tracks
-
-
- fInitialized = kTRUE; //change flag to skip this function next time when processing UserExec
-
-
- fnRCones = TMath::Nint(fRandConeRadiusSquared/fRConesRSquared); //the number of small R=0.1 random cones
-
- // Check for track array
- if(strcmp(fTrackArrayName->Data(), "") != 0){
- fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName->Data()));
- if(!fTrackArray){
- AliWarning(Form("%s: Could not retrieve tracks %s!", GetName(), fTrackArrayName->Data()));
- }else{
- TClass *cl = fTrackArray->GetClass();
- if(!cl->GetBaseClass("AliVParticle")){
- AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data()));
- fTrackArray = 0;
- }
- }
- }
-
- // Check for jet array
- if(strcmp(fJetArrayName->Data(), "") != 0){
- fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName->Data()));
-
- if(!fJetArray){
- AliWarning(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrayName->Data()));
- }else{
- if(!fJetArray->GetClass()->GetBaseClass("AliEmcalJet")){
- AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data()));
- fJetArray = 0;
- }
- }
- }
-
- // Check for list of jets to be removed from background
- if(strcmp(fBackgroundJetArrayName->Data(), "") != 0){
- fBackgroundJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fBackgroundJetArrayName->Data()));
- if(!fBackgroundJetArray){
- AliInfo(Form("%s: Could not retrieve background jets %s!", GetName(), fBackgroundJetArrayName->Data()));
- }else{
- if(!fBackgroundJetArray->GetClass()->GetBaseClass("AliEmcalJet")){
- AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fBackgroundJetArrayName->Data()));
- fBackgroundJetArray = 0;
- }
- }
- }
-
- // Look, if initialization is OK
-
- // Initialize helper class (for vertex selection & pile up correction)
- fHelperClass = new AliAnalysisUtils();
- fHelperClass->SetCutOnZVertexSPD(kFALSE); // kFALSE: no cut; kTRUE: |zvtx-SPD - zvtx-TPC|<0.5cm
-
- return;
-}
-
-
-//________________________________________________________________________
-void AliAnalysisTaskHJetSpectra::GetDeltaPt(Double_t rho1, Double_t &dpt1, Double_t rho2, Double_t &dpt2,
- Double_t rho3, Double_t &dpt3,
- Double_t &rcPhi, Double_t &rcEta,
- Double_t leadingJetExclusionProbability){
-
- //delta pt = random cone - rho
-
- // Define an invalid delta pt
- dpt1 = -10000.0;
- dpt2 = -10000.0;
- dpt3 = -10000.0;
-
- // Define eta range
- Double_t etaMin, etaMax;
- etaMin = -(fTrackEtaWindow-fRandConeRadius);
- etaMax = +(fTrackEtaWindow-fRandConeRadius);
-
- // Define random cone Eta+Phi
- Bool_t coneValid = kTRUE;
- Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
- Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
-
- // if there is a jet, check for overlap if demanded
- if(leadingJetExclusionProbability){
- AliEmcalJet* tmpLeading = NULL;
- Double_t lpt = -1.0;
- // Get leading jet (regardless of pT)
- for(Int_t i = 0; i<fJetArray->GetEntries(); i++){
- AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
- if(!tmpJet) continue;
- if((TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow) && (tmpJet->Area() >= fMinJetArea)){
- if(tmpJet->Pt() > lpt){
- tmpLeading = tmpJet;
- lpt = tmpJet->Pt();
- }
- }
- }
- if(tmpLeading){
- Double_t excludedJetPhi = tmpLeading->Phi();
- Double_t tmpDeltaPhi = RelativePhi(tmpRandConePhi, excludedJetPhi);
- Double_t excludedJetEta = tmpLeading->Eta()-tmpRandConeEta;
-
- // Check, if cone has overlap with jet
- if(tmpDeltaPhi*tmpDeltaPhi + excludedJetEta*excludedJetEta <= fRandConeRadiusSquared){
- // Define probability to exclude the RC
- Double_t probability = leadingJetExclusionProbability;
-
- // Only exclude cone with a given probability
- if(fRandom->Rndm()<=probability) coneValid = kFALSE;
- }
- }
- }
-
- rcPhi = 9999.0;
- rcEta = 9999.0;
-
- // Get the cones' pt and calculate delta pt
- if(coneValid){
- rcPhi = tmpRandConePhi;
- rcEta = tmpRandConeEta;
- Double_t conePt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius);
- dpt1 = conePt - (rho1*fRandConeRadiusSquared*TMath::Pi());
- dpt2 = conePt - (rho2*fRandConeRadiusSquared*TMath::Pi());
- dpt3 = conePt - (rho3*fRandConeRadiusSquared*TMath::Pi());
- }
-
-}
-
-
-//________________________________________________________________________
-void AliAnalysisTaskHJetSpectra::Calculate(AliVEvent* event){
- //Analyze the event and Fill histograms
-
- if(fAnalyzePythia){
- fh1PtHard->Fill(GetPtHard());
- }
-
- if(fAnalyzeHijing){
- fImpParam = GetImpactParameter();
- fhImpactParameter->Fill(fImpParam);
- }
- //_________________________________________________________________
- // FILL EVENT STATISTICS
- fHistEvtSelection->Fill(1); //Count input event
-
- if(!IsEventInAcceptance(event)) return; //post data is in UserExec
-
-
- // Get centrality
- AliCentrality* tmpCentrality = event->GetCentrality();
- if(!tmpCentrality){
- fHistEvtSelection->Fill(4);
- return; //post data is in UserExec
- }
- Double_t centralityPercentile = -1.0;
- Double_t centralityPercentileV0A = -1.0;
- Double_t centralityPercentileV0C = -1.0;
- Double_t centralityPercentileV0M = -1.0;
- Double_t centralityPercentileZNA = -1.0;
- if(tmpCentrality != NULL){
- centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
- centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
- centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
- centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
- centralityPercentileZNA = tmpCentrality->GetCentralityPercentile("ZNA");
- }
-
- if((centralityPercentile < 0.0) || (centralityPercentile > 100.0)){
- AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile));
- fHistEvtSelection->Fill(4);
- return;
- }
- fhCentrality->Fill(centralityPercentile);
- fhCentralityV0M->Fill(centralityPercentileV0M);
- fhCentralityV0A->Fill(centralityPercentileV0A);
- fhCentralityV0C->Fill(centralityPercentileV0C);
- fhCentralityZNA->Fill(centralityPercentileZNA);
-
- fHistEvtSelection->Fill(0); //Count input event
-
- // END EVENT SELECTION
- //___________________________________________________________
-
- //LOOP OVER TRACKS SEARCH FOR TRIGGER
- std::vector<Int_t> trigTracks; //list pf trigger particle indices
- //Bool_t bContainesHighPtTrack = kFALSE;
-
- Int_t nTracks = fTrackArray->GetEntries();
-
- for(Int_t i = 0; i < nTracks; i++){
- AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
-
- if(!track) continue;
-
- if(IsTrackInAcceptance(track)){
- fhTrackPhi->Fill(track->Pt(), RelativePhi(track->Phi(),0.0)); // phi = -pi az pi
- fhTrackEta->Fill(track->Pt(), track->Eta());
- fhTrackCentVsPt->Fill(track->Pt(), centralityPercentile);
-
- if(fTTlow <= track->Pt() && track->Pt() < fTThigh){
- trigTracks.push_back(i); //trigger candidates
- }
-
- //if(track->Pt()>=8.0) bContainesHighPtTrack = kTRUE;
- }
- }
-
- Int_t ntriggers = (Int_t) trigTracks.size();
- Int_t indexSingleRndTrig = -1; //index of single random trigger
- Double_t areaJet, pTJet;
- Double_t tmpArray[3];
- Double_t rhoFromCellMedian = 0.0; //UE density cell median
- Double_t rhoCone = 0.0; //UE density perp cone
- Double_t rhoCMS = 0.0; //UE density ala CMS
- Double_t deltaptCellMedian, deltaptCone, deltaptCMS, randConePhi, randConeEta;
- Double_t distanceFromTrigger;
-
- if(ntriggers>0){
- if(fTTtype==0){ //select single inclusive trigger
- indexSingleRndTrig = fRandom->Integer(ntriggers); //Integer 0 ... ntriggers-1
- }
- }
-
- rhoFromCellMedian = EstimateBgRhoMedian();
- rhoCone = EstimateBgCone();
- rhoCMS = GetExternalRho();
-
- fhRhoCellMedianIncl->Fill((Float_t) rhoFromCellMedian,(Float_t) centralityPercentile);
- fhRhoConeIncl->Fill( (Float_t) rhoCone, (Float_t) centralityPercentile);
- fhRhoCMSIncl->Fill( (Float_t) rhoCMS, (Float_t) centralityPercentile);
-
- for(Int_t irc=0; irc<fNofRandomCones; irc++){ //generate 4 random cones per event
- GetDeltaPt(rhoFromCellMedian, deltaptCellMedian,rhoCone, deltaptCone, rhoCMS, deltaptCMS, randConePhi, randConeEta, 0);
-
- fhDeltaPtMedianIncl->Fill(deltaptCellMedian, (Double_t) centralityPercentile);
- fhDeltaPtConeIncl->Fill( deltaptCone, (Double_t) centralityPercentile);
- fhDeltaPtCMSIncl->Fill( deltaptCMS, (Double_t) centralityPercentile);
-
- if(ntriggers>0){
- //fill delta pt histograms near side + away side
- fhDeltaPtMedian->Fill( deltaptCellMedian, (Double_t) centralityPercentile);
- fhDeltaPtCone->Fill( deltaptCone, (Double_t) centralityPercentile);
- fhDeltaPtCMS->Fill( deltaptCMS, (Double_t) centralityPercentile);
-
- if(indexSingleRndTrig>-1){
- AliVTrack* triggHad = static_cast<AliVTrack*>(fTrackArray->At(trigTracks[indexSingleRndTrig]));
- Double_t dphiTrigRC = RelativePhi(triggHad->Phi(), randConePhi);
- Double_t detaTrigRC = triggHad->Eta()- randConeEta;
- if(TMath::Abs(dphiTrigRC)< TMath::Pi()/2){ //near side
- fhDeltaPtMedianNearSide->Fill( deltaptCellMedian, (Double_t) centralityPercentile);
- fhDeltaPtCMSNearSide->Fill( deltaptCMS, (Double_t) centralityPercentile);
- }else{ //away side
- fhDeltaPtMedianAwaySide->Fill( deltaptCellMedian, (Double_t) centralityPercentile);
- fhDeltaPtCMSAwaySide->Fill( deltaptCMS, (Double_t) centralityPercentile);
- }
-
- distanceFromTrigger = sqrt(dphiTrigRC*dphiTrigRC+detaTrigRC*detaTrigRC);
- while(distanceFromTrigger<0.5 + fRandConeRadius){
- GetDeltaPt(rhoFromCellMedian, deltaptCellMedian,rhoCone, deltaptCone, rhoCMS, deltaptCMS, randConePhi, randConeEta, 0);
- dphiTrigRC = RelativePhi(triggHad->Phi(), randConePhi);
- detaTrigRC = triggHad->Eta()- randConeEta;
- distanceFromTrigger = sqrt(dphiTrigRC*dphiTrigRC+detaTrigRC*detaTrigRC);
- }
- if(distanceFromTrigger>0.5 + fRandConeRadius){
- fhDeltaPtMedianExclTrigCone->Fill( deltaptCellMedian, (Double_t) centralityPercentile);
- fhDeltaPtCMSExclTrigCone->Fill( deltaptCMS, (Double_t) centralityPercentile);
- }
- }
- }
- }
-
- //_______________________________________
- Int_t idxLeadingJetAwaySide = -1;
- Double_t ptLeadingJetAwaySide = -1.0;
- Double_t phiTrigger = -1.0; //-pi,pi
-
- if(ntriggers>0){
- //Estimate UE density
- //Fill once per event
- fhRhoCellMedian->Fill((Float_t) rhoFromCellMedian,(Float_t) centralityPercentile);
- fhRhoCone->Fill( (Float_t) rhoCone, (Float_t) centralityPercentile);
- fhRhoCMS->Fill( (Float_t) rhoCMS, (Float_t) centralityPercentile);
-
-
- //TRIGGER PARTICLE LOOP
- for(Int_t it=0; it<ntriggers; it++){ //loop over trigger configurations
-
- if(fTTtype==0){
- if(it != indexSingleRndTrig) continue;
- }
-
- AliVTrack* triggerHadron = static_cast<AliVTrack*>(fTrackArray->At(trigTracks[it]));
- if(!triggerHadron) continue;
-
- fh2Ntriggers->Fill((Float_t) centralityPercentile, (Float_t) triggerHadron->Pt()); //trigger p
-
- if(fAnalyzeHijing){ //impact parameter for triggered events
- fhImpactParameterTT->Fill(fImpParam);
- }
-
- phiTrigger = RelativePhi(triggerHadron->Phi(),0.0); //-pi,pi
-
- //JET LOOP
- for(Int_t ij = 0; ij < fJetArray->GetEntries(); ij++){
- AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(ij));
- if(!jet){
- AliError(Form("%s: Could not receive jet %d", GetName(), ij));
- continue;
- }
- if(!IsSignalJetInAcceptance(jet)) continue;
-
- areaJet = jet->Area();
- pTJet = jet->Pt();
-
- if(it==0 || it == indexSingleRndTrig){
- fhJetPhi->Fill( pTJet, RelativePhi(jet->Phi(),0.0));
- fhJetEta->Fill( pTJet, jet->Eta());
- }
-
- Double_t dphi = RelativePhi(triggerHadron->Phi(), jet->Phi());
-
- Double_t dfi = dphi; //-0.5*pi to 1.5*Pi
- if(dfi<-0.5*TMath::Pi()) dfi += 2*TMath::Pi();
- if(dfi> 1.5*TMath::Pi()) dfi -= 2*TMath::Pi();
- fhDphiTriggerJetMinBias->Fill((Float_t) jet->Pt(),(Float_t) dfi);
- if(centralityPercentile<20.) fhDphiTriggerJetCent20->Fill((Float_t) jet->Pt(),(Float_t) dfi);
- //-------------------------
-
- if(TMath::Abs(dphi) < fDphiCut) continue; //Dphi cut between trigger and assoc
- fhDphiTriggerJetAccept->Fill(dfi); //Accepted
-
- if(pTJet > ptLeadingJetAwaySide){ //search for the leading away side jet
- idxLeadingJetAwaySide = ij;
- ptLeadingJetAwaySide = pTJet;
- }
-
- //Centrality, A, pTjet
- tmpArray[0] = centralityPercentile;
- tmpArray[1] = areaJet;
- tmpArray[2] = pTJet;
- fHJetSpec->Fill(tmpArray);
-
- //Subtract cell median
- tmpArray[2] = pTJet - areaJet*rhoFromCellMedian;
- fHJetSpecSubUeMedian->Fill(tmpArray);
- fARhoCellMedian->Fill((Float_t) (areaJet*rhoFromCellMedian));
-
- //Subtract perp cone
- tmpArray[2] = pTJet - areaJet*rhoCone;
- fHJetSpecSubUeCone->Fill(tmpArray);
- fARhoCone->Fill((Float_t) (areaJet*rhoCone));
-
- //Subtract CMS bg
- tmpArray[2] = pTJet - areaJet*rhoCMS;
- fHJetSpecSubUeCMS->Fill(tmpArray);
- fARhoCMS->Fill((Float_t) (areaJet*rhoCMS));
-
- }//JET LOOP
- }
- }
-
- //_______________________________________
- // Get delta phi from small R=0.1 cones
- if(ntriggers>0 && indexSingleRndTrig>-1){
-
- AliEmcalJet* jet = NULL;
- Double_t phiExclJet =0., etaExclJet = 9999., rExclJet = 0.0;
-
- if(idxLeadingJetAwaySide>-1){
- jet = static_cast<AliEmcalJet*>(fJetArray->At(idxLeadingJetAwaySide));
- if(!jet){
- AliError(Form("%s: Could not receive leading jet %d", GetName(), idxLeadingJetAwaySide));
- }else{
- phiExclJet = jet->Phi();
- etaExclJet = jet->Eta();
- rExclJet = TMath::Sqrt(jet->Area()/TMath::Pi());
- }
- }
-
- Int_t countPlacedRcones = 0;
- Int_t inwhile=0;
- Double_t dphiMaxFromPiForRcones = TMath::Pi() - fDphiCut + fRandConeRadius - fRConesR;
- Double_t detaMaxForRcones = fTrackEtaWindow - fRConesR;
- Double_t rcphi,rceta;
- Bool_t goodrc;
-
- while(countPlacedRcones < fnRCones){ //generate fnRCones cones of radius R=0.1
- inwhile++;
- if(inwhile>500){
- AliError(Form("%s: Small space where to put another random cone %d", GetName(), idxLeadingJetAwaySide));
- break;
- }
- rcphi = RelativePhi(phiTrigger + TMath::Pi() + dphiMaxFromPiForRcones*fRandom->Uniform(-1.0,1.0), 0.0); //-pi,pi
- rceta = detaMaxForRcones*fRandom->Uniform(-1.0,1.0);
- if(jet){//do not merge random cones with the leading jet
- if(!DistantCones(rcphi,rceta,fRConesR, phiExclJet, etaExclJet, rExclJet)) continue;
- }
-
- goodrc = kTRUE; //generate disjoint random cones
- for(int k=0; k<countPlacedRcones; k++){
- if(!DistantCones(rcphi, rceta, fRConesR, (Double_t) fRConePhi[k], (Double_t) fRConeEta[k], fRConesR)){
- goodrc = kFALSE;
- break;
- }
- }//end loop over already placed cones
-
- if(!goodrc) continue;
- fRConePhi[countPlacedRcones] = rcphi;
- fRConeEta[countPlacedRcones] = rceta;
- countPlacedRcones++;
- }//end of loop generating small R=0.1 random cones
- //sum track pT in the random cones
- Double_t sumPtofRandomCones = 0.0;
-
- for(Int_t itrk = 0; itrk < nTracks; itrk++){
- AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(itrk));
- if(!track) continue;
-
- if(IsTrackInAcceptance(track)){
- for(int k=0; k<countPlacedRcones; k++){
-
- rcphi = RelativePhi(track->Phi(), fRConePhi[k]); // phi = -pi az pi
- rceta = track->Eta() - fRConeEta[k];
- if(rcphi*rcphi + rceta*rceta < fRConesRSquared){
- sumPtofRandomCones += track->Pt(); //track is in the cone
- break;
- }
- }//loop over cones
- }
- }//loop over tracks
- Double_t totarea = countPlacedRcones*TMath::Pi()*fRConesRSquared;
- fhDeltaPtMedianExclAwayJet->Fill( sumPtofRandomCones - rhoFromCellMedian*totarea, (Double_t) centralityPercentile );
- fhDeltaPtCMSExclAwayJet->Fill( sumPtofRandomCones - rhoCMS*totarea , (Double_t) centralityPercentile );
-
- }//end delta phi from small R=0.1 random cones
-
- return;
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskHJetSpectra::UserNotify(){
- // Implemented Notify() to read the cross sections
- // and number of trials from pyxsec.root
- /*
- if(fAnalyzePythia){
-
- TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
- TFile *currFile = tree->GetCurrentFile();
-
- TString file(currFile->GetName());
-
- if(file.Contains("root_archive.zip#")){
- Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
- Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
- file.Replace(pos+1,20,"");
- }
- else {
- // not an archive take the basename....
- file.ReplaceAll(gSystem->BaseName(file.Data()),"");
- }
-
- TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
- if(!fxsec){
- // next trial fetch the histgram file
- fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
- if(!fxsec){
- // not a severe condition but inciate that we have no information
- return kFALSE;
- }
- else{
- // find the tlist we want to be independtent of the name so use the Tkey
- TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
- if(!key){
- fxsec->Close();
- return kFALSE;
- }
- TList *list = dynamic_cast<TList*>(key->ReadObj());
- if(!list){
- fxsec->Close();
- return kFALSE;
- }
- fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
- fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
- fxsec->Close();
- }
- } // no tree pyxsec.root
- else {
- TTree *xtree = (TTree*)fxsec->Get("Xsection");
- if(!xtree){
- fxsec->Close();
- return kFALSE;
- }
- UInt_t ntrials = 0;
- Double_t xsection = 0;
- xtree->SetBranchAddress("xsection",&xsection);
- xtree->SetBranchAddress("ntrials",&ntrials);
- xtree->GetEntry(0);
- fTrials = ntrials;
- fCrossSection = xsection;
- fxsec->Close();
- }
-
-
- fh1Xsec->Fill("<#sigma>", fCrossSection);
- fh1Trials->Fill("#sum{ntrials}",fTrials);
-
- }
- */
- return kTRUE;
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskHJetSpectra::Terminate(Option_t *){
- //Treminate
- PostData(1, fOutputList);
-
- // Mandatory
- fOutputList = dynamic_cast<TList*> (GetOutputData(1)); // '1' refers to the output slot
- if(!fOutputList) {
- printf("ERROR: Output list not available\n");
- return;
- }
-}
-
-//________________________________________________________________________
-AliAnalysisTaskHJetSpectra::~AliAnalysisTaskHJetSpectra(){
- // Destructor. Clean-up the output list, but not the histograms that are put inside
- // (the list is owner and will clean-up these histograms). Protect in PROOF case.
- if(fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
- delete fOutputList;
- }
- delete fRandom;
- delete fTrackArrayName;
- delete fJetArrayName;
- delete fBackgroundJetArrayName;
- delete fHelperClass;
-
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskHJetSpectra::UserCreateOutputObjects(){
- // called once to create user defined output objects like histograms, plots etc.
- // and to put it on the output list.
- // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
-
- fRandom = new TRandom3(0);
-
- fOutputList = new TList();
- fOutputList->SetOwner(); // otherwise it produces leaks in merging
- Bool_t oldStatus = TH1::AddDirectoryStatus();
- TH1::AddDirectory(kFALSE);
-
- //__________________________________________________________
- // Event statistics
- fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
- fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
- fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
- fHistEvtSelection->GetXaxis()->SetBinLabel(3,"pile up (rejected)");
- fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
- fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
-
- fOutputList->Add(fHistEvtSelection);
- //___________________________________________________________
- // Hard trigger counter
- fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
- fNumberOfCentralityBins,0.0,100.0,50,0.0,50.0);
- fOutputList->Add(fh2Ntriggers);
- //___________________________________________________________
- // trigger associated jet spectra (jet pT not corrected for UE)
- Int_t bw = (fUseDoubleBinPrecision==0) ? 1 : 2; //make larger bin width
-
- //jet associated to given TT
- //Centrality, A, pTjet
- const Int_t dimSpec = 3;
- const Int_t nBinsSpec[dimSpec] = {fNumberOfCentralityBins, 50, bw*110};
- const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, -20.0};
- const Double_t hiBinSpec[dimSpec] = {100.0, 1.5, 200.0};
- fHJetSpec = new THnSparseF("fHJetSpec",
- "Recoil jet spectrum [cent,A,pTjet]",
- dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
- fOutputList->Add(fHJetSpec);
- //___________________________________________________________
- //jet associated to given TT (jet pT corrected with rho from cell median)
- fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
- fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]");
- fOutputList->Add(fHJetSpecSubUeMedian);
- //___________________________________________________________
- //jet associated to given TT (jet pT corrected with rho from perp cone)
- fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
- fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]");
- fOutputList->Add(fHJetSpecSubUeCone);
- //___________________________________________________________
- //jet associated to given TT (jet pT corrected with rho from CMS approach)
- fHJetSpecSubUeCMS = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCMS");
- fHJetSpecSubUeCMS->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]");
- fOutputList->Add(fHJetSpecSubUeCMS);
-
- //____________________________________________________________________
- //UE from cell median [Centrality, rho, pTUe ]
-
- fhRhoCellMedian = new TH2F("fhRhoCellMedian","Rho",40, 0.0, 20.0, fNumberOfCentralityBins, 0.0, 100.);
- fOutputList->Add(fhRhoCellMedian);
-
- fhRhoCone = (TH2F*) fhRhoCellMedian->Clone("fhRhoCone");
- fOutputList->Add(fhRhoCone);
-
- fhRhoCMS = (TH2F*) fhRhoCellMedian->Clone("fhRhoCMS");
- fOutputList->Add(fhRhoCMS);
-
- fhRhoCellMedianIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoCellMedianIncl");
- fOutputList->Add(fhRhoCellMedianIncl);
-
- fhRhoConeIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoConeIncl");
- fOutputList->Add(fhRhoConeIncl);
-
- fhRhoCMSIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoCMSIncl");
- fOutputList->Add(fhRhoCMSIncl);
-
- //_______________________________________________________________________
- // rho times area
- fARhoCellMedian = new TH1F("fARhoCellMedian","Area times rho",40, 0.0, 20.0);
- fOutputList->Add(fARhoCellMedian);
-
- fARhoCone = (TH1F*) fARhoCellMedian->Clone("fARhoCone");
- fOutputList->Add(fARhoCone);
-
- fARhoCMS = (TH1F*) fARhoCellMedian->Clone("fARhoCMS");
- fOutputList->Add(fARhoCMS);
-
- //_______________________________________________________________________
- // Delta pt distributions
- fhDeltaPtMedian = new TH2D("fhDeltaPtMedian","DeltaPt", bw*110, -20, 200, fNumberOfCentralityBins,0.0,100.0);
- fOutputList->Add(fhDeltaPtMedian);
-
- fhDeltaPtCone = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCone");
- fOutputList->Add(fhDeltaPtCone);
-
- fhDeltaPtCMS = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMS");
- fOutputList->Add(fhDeltaPtCMS);
-
- fhDeltaPtMedianIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianIncl");
- fOutputList->Add(fhDeltaPtMedianIncl);
-
- fhDeltaPtConeIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtConeIncl");
- fOutputList->Add(fhDeltaPtConeIncl);
-
- fhDeltaPtCMSIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSIncl");
- fOutputList->Add(fhDeltaPtCMSIncl);
-
- fhDeltaPtMedianNearSide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianNearSide");
- fOutputList->Add(fhDeltaPtMedianNearSide);
-
- fhDeltaPtMedianAwaySide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianAwaySide");
- fOutputList->Add(fhDeltaPtMedianAwaySide);
-
- fhDeltaPtCMSNearSide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSNearSide");
- fOutputList->Add(fhDeltaPtCMSNearSide);
-
- fhDeltaPtCMSAwaySide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSAwaySide");
- fOutputList->Add(fhDeltaPtCMSAwaySide);
-
- fhDeltaPtMedianExclTrigCone= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianExclTrigCone");
- fOutputList->Add(fhDeltaPtMedianExclTrigCone);
-
- fhDeltaPtCMSExclTrigCone= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSExclTrigCone");
- fOutputList->Add(fhDeltaPtCMSExclTrigCone);
-
- fhDeltaPtMedianExclAwayJet = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianExclAwayJet");
- fOutputList->Add(fhDeltaPtMedianExclAwayJet);
-
- fhDeltaPtCMSExclAwayJet = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSExclAwayJet");
- fOutputList->Add(fhDeltaPtCMSExclAwayJet);
-
- //_______________________________________________________________________
- //inclusive azimuthal and pseudorapidity histograms
- fhJetPhi = new TH2F("fhJetPhi","Azim dist jets vs pTjet", 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
- fOutputList->Add(fhJetPhi);
- //-------------------------
- fhTrackPhi = new TH2F("fhTrackPhi","azim dist trig had vs pT,trk", 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
- fOutputList->Add(fhTrackPhi);
- //-------------------------
- fhJetEta = new TH2F("fhJetEta","Eta dist jets vs pTjet", 50,0, 100, 40,-0.9,0.9);
- fOutputList->Add(fhJetEta);
- //-------------------------
- fhTrackEta = new TH2F("fhTrackEta","Eta dist trig had vs pT,trk", 50, 0, 50, 40,-0.9,0.9);
- fOutputList->Add(fhTrackEta);
- //-------------------------
- fhTrackCentVsPt = new TH2F("fhTrackCentVsPt","pT,trk vs centrality", 50, 0, 50, fNumberOfCentralityBins,0,100);
- fOutputList->Add(fhTrackCentVsPt);
- //-------------------------
- fhVertexZ = new TH1F("fhVertexZ","z vertex",40,-20,20);
- fOutputList->Add(fhVertexZ);
- //-------------------------
- fhVertexZAccept = new TH1F("fhVertexZAccept","z vertex after cut",40,-20,20);
- fOutputList->Add(fhVertexZAccept);
- //-------------------------
- //fhContribVtx = new TH1F("fhContribVtx","contrib to vtx",200,0,200);
- //fOutputList->Add(fhContribVtx);
- //-------------------------
- //fhContribVtxAccept = new TH1F("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
- //fOutputList->Add(fhContribVtxAccept);
- //-------------------------
- fhDphiTriggerJetMinBias = new TH2F("fhDphiTriggerJetMinBias","Deltaphi trig-jet",50,0,100, 100, -0.5*TMath::Pi(),1.5*TMath::Pi());
- fOutputList->Add(fhDphiTriggerJetMinBias);
-
- fhDphiTriggerJetCent20 = (TH2F*) fhDphiTriggerJetMinBias->Clone("fhDphiTriggerJetCent20");
- fOutputList->Add(fhDphiTriggerJetCent20);
- //-------------------------
-
- fhDphiTriggerJetAccept = new TH1F("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -0.5*TMath::Pi(),1.5*TMath::Pi());
- fOutputList->Add(fhDphiTriggerJetAccept);
- //-------------------------
- fhCentrality = new TH1F("fhCentrality","Centrality",100,0,100);
- fOutputList->Add(fhCentrality);
- //-------------------------
- fhCentralityV0M = new TH1F("hCentralityV0M","hCentralityV0M",100,0,100);
- fOutputList->Add(fhCentralityV0M);
- //-------------------------
- fhCentralityV0A = new TH1F("hCentralityV0A","hCentralityV0A",100,0,100);
- fOutputList->Add(fhCentralityV0A);
- //-------------------------
- fhCentralityV0C = new TH1F("hCentralityV0C","hCentralityV0C",100,0,100);
- fOutputList->Add(fhCentralityV0C);
- //-------------------------
- fhCentralityZNA = new TH1F("hCentralityZNA","hCentralityZNA",100,0,100);
- fOutputList->Add(fhCentralityZNA);
-
- fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
- fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
- fOutputList->Add(fh1Xsec);
-
- fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
- fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
- fOutputList->Add(fh1Trials);
-
- fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",300,0,300);
- fOutputList->Add(fh1PtHard);
-
- fhImpactParameter = new TH1D("fhImpactParameter","impact parameter distribution from HIJING",50,0,10);
- fOutputList->Add(fhImpactParameter);
-
- fhImpactParameterTT = new TH1D("fhImpactParameterTT","b versus TT",50,0,10);
- fOutputList->Add(fhImpactParameterTT);
- // =========== Switch on Sumw2 for all histos ===========
- for(Int_t i=0; i<fOutputList->GetEntries(); i++){
- TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
- if(h1){
- h1->Sumw2();
- continue;
- }
- THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
- if(hn){
- hn->Sumw2();
- }
- }
- TH1::AddDirectory(oldStatus);
-
-
- PostData(1, fOutputList);
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskHJetSpectra::UserExec(Option_t *){
- //executed in each event
-
- if(!InputEvent()){
- AliError("??? Event pointer == 0 ???");
- return;
- }
-
- //Execute only once: Get tracks, jets, background from arrays if not already given
- if(!fInitialized) ExecOnce();
- if(fJetArray && fTrackArray && fBackgroundJetArray){
- Calculate(InputEvent());
- }
- PostData(1, fOutputList);
-}
-//________________________________________________________________________
-
-Double_t AliAnalysisTaskHJetSpectra::RelativePhi(Double_t mphi,Double_t vphi){
- //Get relative azimuthal angle of two particles -pi to pi
- if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
- else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
-
- if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
- else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
-
- Double_t dphi = mphi - vphi;
- if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
- else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
-
- return dphi;//dphi in [-Pi, Pi]
-}
-
-//________________________________________________________________________
-Double_t AliAnalysisTaskHJetSpectra::EstimateBgRhoMedian(){
- //Estimate background rho by means of integrating track pT outside identified jet cones
- Double_t rhoMedian = 0.0;
-
- //phi,eta and R2 of jets to be removed
- std::vector<Double_t> jphi;
- std::vector<Double_t> jeta;
- std::vector<Double_t> jRsquared;
-
- if(!fBackgroundJetArray) return 0.0;
- if(!fTrackArray) return 0.0;
-
- for(Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++){
- AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
-
- if(!backgroundJet){
- AliError(Form("%s: Could not receive jet %d", GetName(), i));
- continue;
- }
- if(!IsBackgroundJetInAcceptance(backgroundJet)) continue; //apply minimum pT cut on jet to be removed from bg
- jphi.push_back(RelativePhi(backgroundJet->Phi(),0.0)); //-pi,pi
- jeta.push_back(backgroundJet->Eta());
- jRsquared.push_back(1.78*backgroundJet->Area()/TMath::Pi()); //1.78 = JetArea_R04/JetArea_R03
- }
-
-
- static Double_t nOutCone[10][4];
- static Double_t sumPtOutOfCone[10][4];
- Double_t rndphi, rndeta;
- Double_t rndphishift, rndetashift;
- Double_t dphi, deta;
- Bool_t bIsInCone;
-
-
- for(Int_t ie=0; ie < fnEta; ie++){
- for(Int_t ip=0; ip < fnPhi; ip++){
- nOutCone[ip][ie] = 0.0; //initialize counter
- sumPtOutOfCone[ip][ie] = 0.0;
- }
- }
-
- //get area in cells out of identified jet cones
- if(jphi.size()==0){ //no jet to be removed from the bg => all areas have their nominal area
- for(Int_t ie=0; ie < fnEta; ie++){
- for(Int_t ip=0; ip < fnPhi; ip++){
- nOutCone[ip][ie] = fNofRndTrials;
- }
- }
- }else{
- for(Int_t it=0; it<fNofRndTrials; it++){
-
- rndphi = fRandom->Uniform(0, fPhiSize);
- rndeta = fRandom->Uniform(0, fEtaSize);
-
- for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
- rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
- for(Int_t ie=0; ie<fnEta; ie++){
- rndetashift = rndeta + ie*fEtaSize - fEtaSize;
-
- bIsInCone = 0; //tag if trial is in the jet cone
- for(Int_t ij=0; ij< (Int_t) jRsquared.size(); ij++){
- deta = jeta[ij] - rndetashift;
- dphi = RelativePhi(rndphishift,jphi[ij]);
- if((dphi*dphi + deta*deta) < jRsquared[ij]){
- bIsInCone = 1;
- break;
- }
- }
- if(!bIsInCone) nOutCone[ip][ie]++;
- }
- }
- }
- }
-
- Int_t phicell,etacell;
- for(Int_t ip=0; ip < fTrackArray->GetEntries(); ip++){
- AliVTrack* part = static_cast<AliVTrack*>(fTrackArray->At(ip));
- if(!part) continue;
- if(!IsTrackInAcceptance((AliVParticle*) part)) continue;
-
- bIsInCone = 0; //init
- for(Int_t ij=0; ij<(Int_t) jRsquared.size(); ij++){
- dphi = RelativePhi(jphi[ij], part->Phi());
- deta = jeta[ij] - part->Eta();
- if((dphi*dphi + deta*deta) < jRsquared[ij]){
- bIsInCone = 1;
- break;
- }
- }
- if(!bIsInCone){
- phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
- etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
- sumPtOutOfCone[phicell][etacell]+= part->Pt();
- }
- }
- // Calculate rho
- static Double_t rhoInCells[20];
- Double_t relativeArea;
- Int_t nCells=0;
- Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
- for(Int_t ip=0; ip<fnPhi; ip++){
- for(Int_t ie=0; ie<fnEta; ie++){
- relativeArea = nOutCone[ip][ie]/fNofRndTrials;
-
- bufferArea += relativeArea;
- bufferPt += sumPtOutOfCone[ip][ie];
- if(bufferArea > fJetFreeAreaFrac){
- rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
-
- bufferArea = 0.0;
- bufferPt = 0.0;
- nCells++;
- }
- }
- }
-
- if(nCells>0){
- rhoMedian = TMath::Median(nCells, rhoInCells);
- }
- return rhoMedian;
-}
-
-//________________________________________________________________________
-Double_t AliAnalysisTaskHJetSpectra::EstimateBgCone(){
- //Estimate background rho by means of integrating track pT outside identified jet cones
- Double_t rhoPerpCone = 0.0;
-
- Double_t pTleading = -1.0;
- Double_t phiLeading = 1000.;
- Double_t etaLeading = 1000.;
-
- if(!fJetArray) return 0.0;
-
- for(Int_t ij = 0; ij < fJetArray->GetEntries(); ij++){
- AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(ij));
- if(!jet){
- AliError(Form("%s: Could not receive jet %d", GetName(), ij));
- continue;
- }
- if(!IsSignalJetInAcceptance(jet)) continue;
-
- if(pTleading < jet->Pt()){
- pTleading = jet->Pt();
- phiLeading = jet->Phi();
- etaLeading = jet->Eta();
- }
- }
- if(pTleading < 0.0) return 0.0;
-
- Double_t phileftcone = phiLeading + TMath::Pi()/2;
- Double_t phirightcone = phiLeading - TMath::Pi()/2;
-
- /* Double_t dp, de;
-
- for(Int_t ip=0; ip < fTrackArray->GetEntries(); ip++){
-
- AliVTrack* part = static_cast<AliVTrack*>(fTrackArray->At(ip));
- if(!part) continue;
- if(!IsTrackInAcceptance((AliVParticle*) part)) continue;
-
-
- dp = RelativePhi(phileftcone, part->Phi());
- de = etaLeading - part->Eta();
- if( dp*dp + de*de < fRandConeRadiusSquared ) rhoPerpCone += part->Pt();
-
- dp = RelativePhi(phirightcone, part->Phi());
- if( dp*dp + de*de < fRandConeRadiusSquared) rhoPerpCone += part->Pt();
-
- }*/
-
- rhoPerpCone += GetConePt(etaLeading, phileftcone, fRandConeRadius);
- rhoPerpCone += GetConePt(etaLeading, phirightcone, fRandConeRadius);
-
- //normalize total pT by two times cone are
- rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fRandConeRadiusSquared);
-
-
-
- return rhoPerpCone;
-}
-//________________________________________________________________________
-Bool_t AliAnalysisTaskHJetSpectra::DistantCones(Double_t phi1, Double_t eta1, Double_t r1, Double_t phi2, Double_t eta2, Double_t r2){
- //checks if the two cones are farther away than the sum of their radii
-
- Double_t dphi = RelativePhi(phi1,phi2);
- Double_t deta = eta1-eta2;
- Double_t d = r1+r2;
- if( dphi*dphi + deta*deta < d*d ) return kFALSE;
-
- return kTRUE;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff