[u/mrichter/AliRoot.git] / PWG / EMCAL / AliAnalysisTaskEmcal.cxx

// $Id: AliAnalysisTaskEmcal.cxx 56756 2012-05-30 05:03:02Z loizides $
//
// Emcal base analysis task.
//
// Author: S.Aiola

#include "AliAnalysisTaskEmcal.h"

#include <TClonesArray.h>
#include <TList.h>
#include <TObject.h>
#include <TH1F.h>

#include "AliAODEvent.h"
#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliEMCALGeometry.h"
#include "AliESDEvent.h"
#include "AliEmcalParticle.h"
#include "AliEventplane.h"
#include "AliInputEventHandler.h"
#include "AliLog.h"
#include "AliMCParticle.h"
#include "AliVCluster.h"
#include "AliVEventHandler.h"
#include "AliVParticle.h"
#include "AliVCaloTrigger.h"

ClassImp(AliAnalysisTaskEmcal)

//________________________________________________________________________
AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() : 
  AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
  fAnaType(kTPC),
  fForceBeamType(kNA),
  fGeneralHistograms(kFALSE),
  fInitialized(kFALSE),
  fCreateHisto(kTRUE),
  fTracksName(),
  fCaloName(),
  fCaloCellsName(),
  fCaloTriggersName(),
  fMinCent(-999),
  fMaxCent(-999),
  fMinVz(-999),
  fMaxVz(-999),
  fOffTrigger(AliVEvent::kAny),
  fTrigClass(),
  fNbins(500),
  fMinBinPt(0),
  fMaxBinPt(250),
  fClusPtCut(0.15),
  fTrackPtCut(0.15),
  fTrackMinEta(-0.9),
  fTrackMaxEta(0.9),
  fTrackMinPhi(-10),
  fTrackMaxPhi(10),
  fClusTimeCutLow(-10),
  fClusTimeCutUp(10),
  fMinPtTrackInEmcal(0),
  fEventPlaneVsEmcal(-1),
  fMinEventPlane(-10),
  fMaxEventPlane(10),
  fCentEst("V0M"),
  fTrackBitMap(0),
  fClusterBitMap(0),
  fMCTrackBitMap(0),
  fMCClusterBitMap(0),
  fMinMCLabel(0),
  fIsEmbedded(kFALSE),
  fNcentBins(4),
  fGeom(0),
  fTracks(0),
  fCaloClusters(0),
  fCaloCells(0),
  fCaloTriggers(0),
  fCent(0),
  fCentBin(-1),
  fEPV0(-1.0),
  fEPV0A(-1.0),
  fEPV0C(-1.0),
  fNVertCont(0),
  fBeamType(kNA),
  fOutput(0),
  fHistCentrality(0),
  fHistZVertex(0),
  fHistEventPlane(0)
{
  // Default constructor.

  fVertex[0] = 0;
  fVertex[1] = 0;
  fVertex[2] = 0;
}

//________________________________________________________________________
AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) : 
  AliAnalysisTaskSE(name),
  fAnaType(kTPC),
  fForceBeamType(kNA),
  fGeneralHistograms(kFALSE),
  fInitialized(kFALSE),
  fCreateHisto(histo),
  fTracksName(),
  fCaloName(),
  fCaloCellsName(),
  fCaloTriggersName(),
  fMinCent(-999),
  fMaxCent(-999),
  fMinVz(-999),
  fMaxVz(-999),
  fOffTrigger(AliVEvent::kAny),
  fTrigClass(),
  fNbins(500),
  fMinBinPt(0),
  fMaxBinPt(250),
  fClusPtCut(0.15),
  fTrackPtCut(0.15),
  fTrackMinEta(-0.9),
  fTrackMaxEta(0.9),
  fTrackMinPhi(-10),
  fTrackMaxPhi(10),
  fClusTimeCutLow(-10),
  fClusTimeCutUp(10),
  fMinPtTrackInEmcal(0),
  fEventPlaneVsEmcal(-1),
  fMinEventPlane(-10),
  fMaxEventPlane(10),
  fCentEst("V0M"),
  fTrackBitMap(0),
  fClusterBitMap(0),
  fMCTrackBitMap(0),
  fMCClusterBitMap(0),
  fMinMCLabel(0),
  fIsEmbedded(kFALSE),
  fNcentBins(4),
  fGeom(0),
  fTracks(0),
  fCaloClusters(0),
  fCaloCells(0),
  fCaloTriggers(0),
  fCent(0),
  fCentBin(-1),
  fEPV0(-1.0),
  fEPV0A(-1.0),
  fEPV0C(-1.0),
  fNVertCont(0),
  fBeamType(kNA),
  fOutput(0),
  fHistCentrality(0),
  fHistZVertex(0),
  fHistEventPlane(0)
{
  // Standard constructor.

  fVertex[0] = 0;
  fVertex[1] = 0;
  fVertex[2] = 0;

  if (fCreateHisto) {
    DefineOutput(1, TList::Class()); 
  }
}

//________________________________________________________________________
AliAnalysisTaskEmcal::~AliAnalysisTaskEmcal()
{
  // Destructor
}

//________________________________________________________________________
void AliAnalysisTaskEmcal::UserCreateOutputObjects()
{
  // Create user output.
  if (!fCreateHisto)
    return;

  OpenFile(1);
  fOutput = new TList();
  fOutput->SetOwner();

  if (fForceBeamType == kpp)
    fNcentBins = 1;

  if (!fGeneralHistograms)
    return;

  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
  fHistCentrality->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistCentrality);

  fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
  fHistZVertex->GetXaxis()->SetTitle("z");
  fHistZVertex->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistZVertex);

  fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
  fHistEventPlane->GetXaxis()->SetTitle("event plane");
  fHistEventPlane->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistEventPlane);

  PostData(1, fOutput);
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcal::FillGeneralHistograms()
{
  fHistCentrality->Fill(fCent);
  fHistZVertex->Fill(fVertex[2]);
  fHistEventPlane->Fill(fEPV0);

  return kTRUE;
}

//________________________________________________________________________
void AliAnalysisTaskEmcal::UserExec(Option_t *) 
{
  // Main loop, called for each event.

  if (!fInitialized)
    ExecOnce();

  if (!fInitialized)
    return;

  if (!RetrieveEventObjects())
    return;

  if (!IsEventSelected()) 
    return;

  if (fGeneralHistograms && fCreateHisto) {
    if (!FillGeneralHistograms())
      return;
  }

  if (!Run())
    return;

  if (fCreateHisto) {
    if (!FillHistograms())
      return;
  }
    
  if (fCreateHisto && fOutput) {
    // information for this iteration of the UserExec in the container
    PostData(1, fOutput);
  }
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus) const
{
  // Return true if cluster is accepted.

  if (!clus)
    return kFALSE;
      
  if (!clus->IsEMCAL())
    return kFALSE;

  if (clus->GetLabel() > fMinMCLabel) {
    if (clus->TestBits(fMCClusterBitMap) != (Int_t)fMCClusterBitMap) {
      AliDebug(2,"MC Cluster not accepted because of bit map.");
      return kFALSE;
    }
  }
  else {
    if (clus->TestBits(fClusterBitMap) != (Int_t)fClusterBitMap) {
      AliDebug(2,"Cluster not accepted because of bit map.");
      return kFALSE;
    }
  }

  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow)
    return kFALSE;

  TLorentzVector nPart;
  clus->GetMomentum(nPart, const_cast<Double_t*>(fVertex));

  if (nPart.Et() < fClusPtCut)
    return kFALSE;

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track) const
{
  // Return true if track is accepted.

  if (!track)
    return kFALSE;

  if (TMath::Abs(track->GetLabel()) > fMinMCLabel) {
    if(track->TestBits(fMCTrackBitMap) != (Int_t)fMCTrackBitMap) {
      AliDebug(2,"MC Track not accepted because of bit map.");
      return kFALSE;
    }
  }
  else {
    if(track->TestBits(fTrackBitMap) != (Int_t)fTrackBitMap) {
      AliDebug(2,"Track not accepted because of bit map.");
      return kFALSE;
    }
  }

  if (track->Pt() < fTrackPtCut)
    return kFALSE;

  if (track->Eta() < fTrackMinEta || track->Eta() > fTrackMaxEta || 
      track->Phi() < fTrackMinPhi || track->Phi() > fTrackMaxPhi)
    return kFALSE;
  
  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcal::AcceptEmcalPart(AliEmcalParticle *part) const
{
  // Return true if EMCal particle is accepted.

  if (!part)
    return kFALSE;

  if (part->IsTrack()) {
    if (part->IsMC()) { 
      if (part->TestBits(fMCTrackBitMap) != (Int_t)fMCTrackBitMap)
	return kFALSE;
    }
    else {
      if (part->TestBits(fTrackBitMap) != (Int_t)fTrackBitMap)
	return kFALSE;
    }

    if (part->Pt() < fTrackPtCut)
      return kFALSE;

    if (part->Eta() < fTrackMinEta || part->Eta() > fTrackMaxEta || 
	part->Phi() < fTrackMinPhi || part->Phi() > fTrackMaxPhi)
      return kFALSE;
  }

  if (part->IsCluster()) {
    if (part->IsMC(fMinMCLabel)) { 
      if (part->TestBits(fMCClusterBitMap) != (Int_t)fMCClusterBitMap)
	return kFALSE;
    }
    else {
      if (part->TestBits(fClusterBitMap) != (Int_t)fClusterBitMap)
	return kFALSE;
    }

    if (!part->IsEMCAL())
      return kFALSE;

    if (part->Pt() < fClusPtCut)
      return kFALSE;
  }

  return kTRUE;
}

//________________________________________________________________________
void AliAnalysisTaskEmcal::ExecOnce()
{
  // Init the analysis.

  if (!InputEvent()) {
    AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
    return;
  }

  fGeom = AliEMCALGeometry::GetInstance();
  if (!fGeom) {
    AliError(Form("%s: Can not create geometry", GetName()));
    return;
  }

  if (fEventPlaneVsEmcal >= 0) {
    Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
    fMinEventPlane = ep - TMath::Pi() / 4;
    fMaxEventPlane = ep + TMath::Pi() / 4;
  }

  if (!fCaloName.IsNull() && !fCaloClusters) {
    fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
    if (!fCaloClusters) {
      AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCaloName.Data())); 
      return;
    } else {
      TClass *cl = fCaloClusters->GetClass();
      if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
	AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCaloName.Data())); 
	fCaloClusters = 0;
	return;
      }
    }
  }

  if (!fTracksName.IsNull() && !fTracks) {
    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
    if (!fTracks) {
      AliError(Form("%s: Could not retrieve tracks %s!", GetName(), fTracksName.Data())); 
      return;
    } else {
      TClass *cl = fTracks->GetClass();
      if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
	AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracksName.Data())); 
	fTracks = 0;
	return;
      }
    }
  }

  if (!fCaloCellsName.IsNull() && !fCaloCells) {
    fCaloCells =  dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
    if (!fCaloCells) {
      AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data())); 
      return;
    }
  }

  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
    fCaloTriggers =  dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
    if (!fCaloTriggers) {
      AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data())); 
      return;
    }
  }

  fInitialized = kTRUE;
}

//_____________________________________________________
AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
{
  // Get beam type : pp-AA-pA
  // ESDs have it directly, AODs get it from hardcoded run number ranges

  if (fForceBeamType != kNA)
    return fForceBeamType;

  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
  if (esd) {
    const AliESDRun *run = esd->GetESDRun();
    TString beamType = run->GetBeamType();
    if (beamType == "p-p")
      return kpp;
    else if (beamType == "A-A")
      return kAA;
    else if (beamType == "p-A")
      return kpA;
    else
      return kNA;
  } else {
    Int_t runNumber = InputEvent()->GetRunNumber();
    if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
	(runNumber >= 166529 && runNumber <= 170593))    // LHC11h
    { 
      return kAA;
    } 
    else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
	     (runNumber >= 195344 && runNumber <= 196608))  // LHC13b-f
    {
      return kpA;
    } else {
      return kpp;
    }
  }  
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcal::IsEventSelected()
{
  // Check if event is selected

  if (fOffTrigger != AliVEvent::kAny) {
    UInt_t res = 0;
    const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
    if (eev) {
      res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
    } else {
      const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
      if (aev) {
        res = aev->GetHeader()->GetOfflineTrigger();
      }
    }
    if ((res & fOffTrigger) == 0)
      return kFALSE;
  }

  if (!fTrigClass.IsNull()) {
    TString fired;
    const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
    if (eev) {
      fired = eev->GetFiredTriggerClasses();
    } else {
      const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
      if (aev) {
        fired = aev->GetFiredTriggerClasses();
      }
    }
    if (!fired.Contains("-B-"))
      return kFALSE;
    TObjArray *arr = fTrigClass.Tokenize("|");
    if (!arr)
      return kFALSE;
    Bool_t match = 0;
    for (Int_t i=0;i<arr->GetEntriesFast();++i) {
      TObject *obj = arr->At(i);
      if (!obj)
        continue;
      if (fired.Contains(obj->GetName())) {
        match = 1;
        break;
      }
    }
    delete arr;
    if (!match)
      return kFALSE;
  }

  if ((fMinCent != -999) && (fMaxCent != -999)) {
    if (fCent<fMinCent)
      return kFALSE;
    if (fCent>fMaxCent)
      return kFALSE;
  }

  if ((fMinVz != -999) && (fMaxVz != -999)) {
    if (fNVertCont == 0 )
      return kFALSE;
    Double_t vz = fVertex[2];
    if (vz<fMinVz)
      return kFALSE;
    if (vz>fMaxVz)
      return kFALSE;
  }

  if (fMinPtTrackInEmcal > 0 && fTracks && fGeom) {
    Bool_t trackInEmcalOk = kFALSE;
    Int_t ntracks = fTracks->GetEntries();
    for (Int_t i = 0; i < ntracks; i++) {
      AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(i));
      if (!AcceptTrack(track))
	continue;
      if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() ||
	  track->Phi() < fGeom->GetArm1PhiMin() * TMath::DegToRad() || track->Phi() > fGeom->GetArm1PhiMax() * TMath::DegToRad())
	continue;
      if (track->Pt() > fMinPtTrackInEmcal) {
	trackInEmcalOk = kTRUE;
	break;
      }
    }
    if (!trackInEmcalOk)
      return kFALSE;
  }


  if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
      !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
      !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane)) 
    return kFALSE;


  return kTRUE;
}

//________________________________________________________________________
TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
{
  // Get array from event.

  TClonesArray *arr = 0;
  TString sname(name);
  if (!sname.IsNull()) {
    arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
    if (!arr) {
      AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name)); 
      return 0;
    }
  } else {
    return 0;
  }

  if (!clname)
    return arr;

  TString objname(arr->GetClass()->GetName());
  TClass cls(objname);
  if (!cls.InheritsFrom(clname)) {
    AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!", 
                    GetName(), cls.GetName(), name, clname)); 
    return 0;
  }
  return arr;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
{
  // Retrieve objects from event.

  fVertex[0] = 0;
  fVertex[1] = 0;
  fVertex[2] = 0;
  fNVertCont = 0;

  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
  if (vert) {
    vert->GetXYZ(fVertex);
    fNVertCont = vert->GetNContributors();
  }

  fBeamType = GetBeamType();

  if (fBeamType == kAA || fBeamType == kpA ) {
    AliCentrality *aliCent = InputEvent()->GetCentrality();
    if (aliCent) {
      fCent = aliCent->GetCentralityPercentile(fCentEst.Data()); 
      if(fNcentBins==4) {
	if      (fCent >=  0 && fCent <   10) fCentBin = 0;
	else if (fCent >= 10 && fCent <   30) fCentBin = 1;
	else if (fCent >= 30 && fCent <   50) fCentBin = 2;
	else if (fCent >= 50 && fCent <= 100) fCentBin = 3; 
	else {
	  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
	  fCentBin = 3;
	}
      }
      else {
	Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
	fCentBin = TMath::FloorNint(fCent/centWidth);
      }
    } else {
      AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
      fCentBin = 3;
    }
    AliEventplane *aliEP = InputEvent()->GetEventplane();
    if (aliEP) {
      fEPV0  = aliEP->GetEventplane("V0" ,InputEvent());
      fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
      fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
    } else {
      AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
    }
  } else {
    fCent = 99;
    fCentBin = 0;
  }

  return kTRUE;
}
Commit	Line	Data
2da09763	1	// $Id: AliAnalysisTaskEmcal.cxx 56756 2012-05-30 05:03:02Z loizides $
	2	//
	3	// Emcal base analysis task.
	4	//
	5	// Author: S.Aiola
	6
	7	#include "AliAnalysisTaskEmcal.h"
	8
2da09763	9	#include <TClonesArray.h>
2da09763	10	#include <TList.h>
09ca054b	11	#include <TObject.h>
159802e2	12	#include <TH1F.h>
2da09763	13
d29ed625	14	#include "AliAODEvent.h"
2da09763	15	#include "AliAnalysisManager.h"
2da09763	16	#include "AliCentrality.h"
09ca054b	17	#include "AliEMCALGeometry.h"
09ca054b	18	#include "AliESDEvent.h"
09ca054b	19	#include "AliEmcalParticle.h"
d29ed625	20	#include "AliEventplane.h"
d29ed625	21	#include "AliInputEventHandler.h"
2da09763	22	#include "AliLog.h"
09ca054b	23	#include "AliMCParticle.h"
	24	#include "AliVCluster.h"
	25	#include "AliVEventHandler.h"
	26	#include "AliVParticle.h"
00c77045	27	#include "AliVCaloTrigger.h"
2da09763	28
	29	ClassImp(AliAnalysisTaskEmcal)
	30
	31	//________________________________________________________________________
	32	AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() :
	33	AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
	34	fAnaType(kTPC),
159802e2	35	fForceBeamType(kNA),
159802e2	36	fGeneralHistograms(kFALSE),
2da09763	37	fInitialized(kFALSE),
2da09763	38	fCreateHisto(kTRUE),
20f79d37	39	fTracksName(),
20f79d37	40	fCaloName(),
f660c2d6	41	fCaloCellsName(),
00c77045	42	fCaloTriggersName(),
2ac8ca4f	43	fMinCent(-999),
	44	fMaxCent(-999),
	45	fMinVz(-999),
	46	fMaxVz(-999),
	47	fOffTrigger(AliVEvent::kAny),
b4339010	48	fTrigClass(),
2da09763	49	fNbins(500),
	50	fMinBinPt(0),
	51	fMaxBinPt(250),
e44e8726	52	fClusPtCut(0.15),
e44e8726	53	fTrackPtCut(0.15),
159802e2	54	fTrackMinEta(-0.9),
	55	fTrackMaxEta(0.9),
	56	fTrackMinPhi(-10),
	57	fTrackMaxPhi(10),
8210e4cb	58	fClusTimeCutLow(-10),
8210e4cb	59	fClusTimeCutUp(10),
6c41fc8c	60	fMinPtTrackInEmcal(0),
6f6dc2ba	61	fEventPlaneVsEmcal(-1),
	62	fMinEventPlane(-10),
	63	fMaxEventPlane(10),
a52de945	64	fCentEst("V0M"),
2103dc6a	65	fTrackBitMap(0),
2103dc6a	66	fClusterBitMap(0),
507f74bc	67	fMCTrackBitMap(0),
507f74bc	68	fMCClusterBitMap(0),
5ded2286	69	fMinMCLabel(0),
8c0d179d	70	fIsEmbedded(kFALSE),
2130abba	71	fNcentBins(4),
159802e2	72	fGeom(0),
2da09763	73	fTracks(0),
2da09763	74	fCaloClusters(0),
f660c2d6	75	fCaloCells(0),
00c77045	76	fCaloTriggers(0),
2da09763	77	fCent(0),
2da09763	78	fCentBin(-1),
3c124212	79	fEPV0(-1.0),
	80	fEPV0A(-1.0),
	81	fEPV0C(-1.0),
b4339010	82	fNVertCont(0),
1f6fff78	83	fBeamType(kNA),
a4e91fed	84	fOutput(0),
a4e91fed	85	fHistCentrality(0),
6f6dc2ba	86	fHistZVertex(0),
6f6dc2ba	87	fHistEventPlane(0)
2da09763	88	{
	89	// Default constructor.
	90
	91	fVertex[0] = 0;
	92	fVertex[1] = 0;
	93	fVertex[2] = 0;
2da09763	94	}
	95
	96	//________________________________________________________________________
	97	AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
	98	AliAnalysisTaskSE(name),
	99	fAnaType(kTPC),
159802e2	100	fForceBeamType(kNA),
159802e2	101	fGeneralHistograms(kFALSE),
2da09763	102	fInitialized(kFALSE),
2da09763	103	fCreateHisto(histo),
20f79d37	104	fTracksName(),
20f79d37	105	fCaloName(),
f660c2d6	106	fCaloCellsName(),
00c77045	107	fCaloTriggersName(),
2ac8ca4f	108	fMinCent(-999),
	109	fMaxCent(-999),
	110	fMinVz(-999),
	111	fMaxVz(-999),
	112	fOffTrigger(AliVEvent::kAny),
b4339010	113	fTrigClass(),
2da09763	114	fNbins(500),
	115	fMinBinPt(0),
	116	fMaxBinPt(250),
e44e8726	117	fClusPtCut(0.15),
e44e8726	118	fTrackPtCut(0.15),
159802e2	119	fTrackMinEta(-0.9),
	120	fTrackMaxEta(0.9),
	121	fTrackMinPhi(-10),
	122	fTrackMaxPhi(10),
8210e4cb	123	fClusTimeCutLow(-10),
8210e4cb	124	fClusTimeCutUp(10),
6c41fc8c	125	fMinPtTrackInEmcal(0),
6f6dc2ba	126	fEventPlaneVsEmcal(-1),
	127	fMinEventPlane(-10),
	128	fMaxEventPlane(10),
a52de945	129	fCentEst("V0M"),
2103dc6a	130	fTrackBitMap(0),
2103dc6a	131	fClusterBitMap(0),
507f74bc	132	fMCTrackBitMap(0),
507f74bc	133	fMCClusterBitMap(0),
5ded2286	134	fMinMCLabel(0),
8c0d179d	135	fIsEmbedded(kFALSE),
2130abba	136	fNcentBins(4),
159802e2	137	fGeom(0),
2da09763	138	fTracks(0),
2da09763	139	fCaloClusters(0),
f660c2d6	140	fCaloCells(0),
00c77045	141	fCaloTriggers(0),
2da09763	142	fCent(0),
2da09763	143	fCentBin(-1),
3c124212	144	fEPV0(-1.0),
	145	fEPV0A(-1.0),
	146	fEPV0C(-1.0),
b4339010	147	fNVertCont(0),
1f6fff78	148	fBeamType(kNA),
a4e91fed	149	fOutput(0),
a4e91fed	150	fHistCentrality(0),
6f6dc2ba	151	fHistZVertex(0),
6f6dc2ba	152	fHistEventPlane(0)
2da09763	153	{
	154	// Standard constructor.
	155
	156	fVertex[0] = 0;
	157	fVertex[1] = 0;
	158	fVertex[2] = 0;
	159
	160	if (fCreateHisto) {
2da09763	161	DefineOutput(1, TList::Class());
	162	}
	163	}
	164
	165	//________________________________________________________________________
	166	AliAnalysisTaskEmcal::~AliAnalysisTaskEmcal()
	167	{
	168	// Destructor
	169	}
	170
159802e2	171	//________________________________________________________________________
	172	void AliAnalysisTaskEmcal::UserCreateOutputObjects()
	173	{
	174	// Create user output.
	175	if (!fCreateHisto)
	176	return;
	177
	178	OpenFile(1);
	179	fOutput = new TList();
	180	fOutput->SetOwner();
	181
2130abba	182	if (fForceBeamType == kpp)
	183	fNcentBins = 1;
	184
159802e2	185	if (!fGeneralHistograms)
	186	return;
	187
	188	fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
	189	fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
	190	fHistCentrality->GetYaxis()->SetTitle("counts");
	191	fOutput->Add(fHistCentrality);
	192
	193	fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
	194	fHistZVertex->GetXaxis()->SetTitle("z");
	195	fHistZVertex->GetYaxis()->SetTitle("counts");
	196	fOutput->Add(fHistZVertex);
	197
6f6dc2ba	198	fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
	199	fHistEventPlane->GetXaxis()->SetTitle("event plane");
	200	fHistEventPlane->GetYaxis()->SetTitle("counts");
	201	fOutput->Add(fHistEventPlane);
	202
159802e2	203	PostData(1, fOutput);
	204	}
	205
	206	//________________________________________________________________________
	207	Bool_t AliAnalysisTaskEmcal::FillGeneralHistograms()
	208	{
	209	fHistCentrality->Fill(fCent);
	210	fHistZVertex->Fill(fVertex[2]);
6f6dc2ba	211	fHistEventPlane->Fill(fEPV0);
159802e2	212
	213	return kTRUE;
	214	}
	215
2da09763	216	//________________________________________________________________________
629c7ac0	217	void AliAnalysisTaskEmcal::UserExec(Option_t *)
	218	{
	219	// Main loop, called for each event.
	220
e44e8726	221	if (!fInitialized)
b3376347	222	ExecOnce();
629c7ac0	223
e44e8726	224	if (!fInitialized)
	225	return;
	226
629c7ac0	227	if (!RetrieveEventObjects())
	228	return;
	229
2ac8ca4f	230	if (!IsEventSelected())
	231	return;
	232
159802e2	233	if (fGeneralHistograms && fCreateHisto) {
	234	if (!FillGeneralHistograms())
	235	return;
	236	}
	237
629c7ac0	238	if (!Run())
	239	return;
	240
159802e2	241	if (fCreateHisto) {
	242	if (!FillHistograms())
	243	return;
	244	}
629c7ac0	245
2ac8ca4f	246	if (fCreateHisto && fOutput) {
629c7ac0	247	// information for this iteration of the UserExec in the container
	248	PostData(1, fOutput);
	249	}
	250	}
	251
	252	//________________________________________________________________________
2103dc6a	253	Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus) const
2da09763	254	{
629c7ac0	255	// Return true if cluster is accepted.
	256
	257	if (!clus)
	258	return kFALSE;
fde82e42	259
	260	if (!clus->IsEMCAL())
	261	return kFALSE;
629c7ac0	262
5ded2286	263	if (clus->GetLabel() > fMinMCLabel) {
507f74bc	264	if (clus->TestBits(fMCClusterBitMap) != (Int_t)fMCClusterBitMap) {
	265	AliDebug(2,"MC Cluster not accepted because of bit map.");
	266	return kFALSE;
	267	}
2103dc6a	268	}
507f74bc	269	else {
	270	if (clus->TestBits(fClusterBitMap) != (Int_t)fClusterBitMap) {
	271	AliDebug(2,"Cluster not accepted because of bit map.");
	272	return kFALSE;
	273	}
	274	}
629c7ac0	275
090a0c3e	276	if (clus->GetTOF() > fClusTimeCutUp \|\| clus->GetTOF() < fClusTimeCutLow)
	277	return kFALSE;
	278
629c7ac0	279	TLorentzVector nPart;
	280	clus->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
	281
e44e8726	282	if (nPart.Et() < fClusPtCut)
629c7ac0	283	return kFALSE;
	284
	285	return kTRUE;
	286	}
	287
	288	//________________________________________________________________________
2103dc6a	289	Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track) const
629c7ac0	290	{
159802e2	291	// Return true if track is accepted.
629c7ac0	292
159802e2	293	if (!track)
629c7ac0	294	return kFALSE;
629c7ac0	295
5ded2286	296	if (TMath::Abs(track->GetLabel()) > fMinMCLabel) {
507f74bc	297	if(track->TestBits(fMCTrackBitMap) != (Int_t)fMCTrackBitMap) {
	298	AliDebug(2,"MC Track not accepted because of bit map.");
	299	return kFALSE;
	300	}
	301	}
	302	else {
	303	if(track->TestBits(fTrackBitMap) != (Int_t)fTrackBitMap) {
	304	AliDebug(2,"Track not accepted because of bit map.");
	305	return kFALSE;
	306	}
2103dc6a	307	}
629c7ac0	308
159802e2	309	if (track->Pt() < fTrackPtCut)
629c7ac0	310	return kFALSE;
629c7ac0	311
159802e2	312	if (track->Eta() < fTrackMinEta \|\| track->Eta() > fTrackMaxEta \|\|
159802e2	313	track->Phi() < fTrackMinPhi \|\| track->Phi() > fTrackMaxPhi)
629c7ac0	314	return kFALSE;
159802e2	315
629c7ac0	316	return kTRUE;
	317	}
	318
	319	//________________________________________________________________________
2103dc6a	320	Bool_t AliAnalysisTaskEmcal::AcceptEmcalPart(AliEmcalParticle *part) const
629c7ac0	321	{
159802e2	322	// Return true if EMCal particle is accepted.
629c7ac0	323
159802e2	324	if (!part)
629c7ac0	325	return kFALSE;
629c7ac0	326
2103dc6a	327	if (part->IsTrack()) {
507f74bc	328	if (part->IsMC()) {
	329	if (part->TestBits(fMCTrackBitMap) != (Int_t)fMCTrackBitMap)
	330	return kFALSE;
	331	}
	332	else {
	333	if (part->TestBits(fTrackBitMap) != (Int_t)fTrackBitMap)
	334	return kFALSE;
	335	}
2103dc6a	336
159802e2	337	if (part->Pt() < fTrackPtCut)
159802e2	338	return kFALSE;
629c7ac0	339
159802e2	340	if (part->Eta() < fTrackMinEta \|\| part->Eta() > fTrackMaxEta \|\|
	341	part->Phi() < fTrackMinPhi \|\| part->Phi() > fTrackMaxPhi)
	342	return kFALSE;
	343	}
	344
	345	if (part->IsCluster()) {
5ded2286	346	if (part->IsMC(fMinMCLabel)) {
507f74bc	347	if (part->TestBits(fMCClusterBitMap) != (Int_t)fMCClusterBitMap)
	348	return kFALSE;
	349	}
	350	else {
	351	if (part->TestBits(fClusterBitMap) != (Int_t)fClusterBitMap)
	352	return kFALSE;
	353	}
2103dc6a	354
159802e2	355	if (!part->IsEMCAL())
159802e2	356	return kFALSE;
9733b37f	357
159802e2	358	if (part->Pt() < fClusPtCut)
	359	return kFALSE;
	360	}
	361
629c7ac0	362	return kTRUE;
2da09763	363	}
2da09763	364
629c7ac0	365	//________________________________________________________________________
b3376347	366	void AliAnalysisTaskEmcal::ExecOnce()
8a28ded1	367	{
8a28ded1	368	// Init the analysis.
2ac8ca4f	369
b3376347	370	if (!InputEvent()) {
	371	AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
	372	return;
	373	}
	374
159802e2	375	fGeom = AliEMCALGeometry::GetInstance();
bb98e333	376	if (!fGeom) {
	377	AliError(Form("%s: Can not create geometry", GetName()));
	378	return;
	379	}
159802e2	380
6f6dc2ba	381	if (fEventPlaneVsEmcal >= 0) {
	382	Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
	383	fMinEventPlane = ep - TMath::Pi() / 4;
	384	fMaxEventPlane = ep + TMath::Pi() / 4;
	385	}
	386
159802e2	387	if (!fCaloName.IsNull() && !fCaloClusters) {
b3376347	388	fCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
	389	if (!fCaloClusters) {
	390	AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCaloName.Data()));
	391	return;
	392	} else {
	393	TClass *cl = fCaloClusters->GetClass();
	394	if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
	395	AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCaloName.Data()));
	396	fCaloClusters = 0;
	397	return;
	398	}
	399	}
	400	}
	401
159802e2	402	if (!fTracksName.IsNull() && !fTracks) {
b3376347	403	fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
	404	if (!fTracks) {
	405	AliError(Form("%s: Could not retrieve tracks %s!", GetName(), fTracksName.Data()));
	406	return;
	407	} else {
	408	TClass *cl = fTracks->GetClass();
	409	if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
	410	AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracksName.Data()));
	411	fTracks = 0;
	412	return;
	413	}
	414	}
	415	}
159802e2	416
f660c2d6	417	if (!fCaloCellsName.IsNull() && !fCaloCells) {
	418	fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
	419	if (!fCaloCells) {
507f74bc	420	AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
f660c2d6	421	return;
	422	}
	423	}
	424
00c77045	425	if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
	426	fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
	427	if (!fCaloTriggers) {
	428	AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
	429	return;
	430	}
	431	}
	432
159802e2	433	fInitialized = kTRUE;
8a28ded1	434	}
8a28ded1	435
2ac8ca4f	436	//_____________________________________________________
	437	AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
	438	{
	439	// Get beam type : pp-AA-pA
	440	// ESDs have it directly, AODs get it from hardcoded run number ranges
	441
159802e2	442	if (fForceBeamType != kNA)
	443	return fForceBeamType;
	444
2ac8ca4f	445	AliESDEvent esd = dynamic_cast<AliESDEvent>(InputEvent());
	446	if (esd) {
	447	const AliESDRun *run = esd->GetESDRun();
	448	TString beamType = run->GetBeamType();
	449	if (beamType == "p-p")
	450	return kpp;
	451	else if (beamType == "A-A")
	452	return kAA;
	453	else if (beamType == "p-A")
	454	return kpA;
	455	else
	456	return kNA;
	457	} else {
	458	Int_t runNumber = InputEvent()->GetRunNumber();
	459	if ((runNumber >= 136851 && runNumber <= 139517) \|\| // LHC10h
	460	(runNumber >= 166529 && runNumber <= 170593)) // LHC11h
493d6e49	461	{
2ac8ca4f	462	return kAA;
493d6e49	463	}
	464	else if ((runNumber>=188365 && runNumber <= 188366) \|\| // LHC12g
	465	(runNumber >= 195344 && runNumber <= 196608)) // LHC13b-f
	466	{
	467	return kpA;
2ac8ca4f	468	} else {
	469	return kpp;
	470	}
	471	}
	472	}
	473
	474	//________________________________________________________________________
d29ed625	475	Bool_t AliAnalysisTaskEmcal::IsEventSelected()
2ac8ca4f	476	{
	477	// Check if event is selected
	478
	479	if (fOffTrigger != AliVEvent::kAny) {
d29ed625	480	UInt_t res = 0;
	481	const AliESDEvent eev = dynamic_cast<const AliESDEvent>(InputEvent());
	482	if (eev) {
	483	res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
	484	} else {
	485	const AliAODEvent aev = dynamic_cast<const AliAODEvent>(InputEvent());
	486	if (aev) {
	487	res = aev->GetHeader()->GetOfflineTrigger();
	488	}
	489	}
b4339010	490	if ((res & fOffTrigger) == 0)
	491	return kFALSE;
	492	}
	493
	494	if (!fTrigClass.IsNull()) {
	495	TString fired;
	496	const AliESDEvent eev = dynamic_cast<const AliESDEvent>(InputEvent());
	497	if (eev) {
	498	fired = eev->GetFiredTriggerClasses();
	499	} else {
	500	const AliAODEvent aev = dynamic_cast<const AliAODEvent>(InputEvent());
	501	if (aev) {
	502	fired = aev->GetFiredTriggerClasses();
	503	}
	504	}
	505	if (!fired.Contains("-B-"))
	506	return kFALSE;
	507	TObjArray *arr = fTrigClass.Tokenize("\|");
	508	if (!arr)
	509	return kFALSE;
	510	Bool_t match = 0;
	511	for (Int_t i=0;i<arr->GetEntriesFast();++i) {
	512	TObject *obj = arr->At(i);
	513	if (!obj)
	514	continue;
	515	if (fired.Contains(obj->GetName())) {
	516	match = 1;
	517	break;
	518	}
	519	}
	520	delete arr;
	521	if (!match)
2ac8ca4f	522	return kFALSE;
	523	}
	524
	525	if ((fMinCent != -999) && (fMaxCent != -999)) {
	526	if (fCent<fMinCent)
	527	return kFALSE;
	528	if (fCent>fMaxCent)
	529	return kFALSE;
	530	}
	531
	532	if ((fMinVz != -999) && (fMaxVz != -999)) {
b4339010	533	if (fNVertCont == 0 )
b4339010	534	return kFALSE;
2ac8ca4f	535	Double_t vz = fVertex[2];
	536	if (vz<fMinVz)
	537	return kFALSE;
	538	if (vz>fMaxVz)
	539	return kFALSE;
	540	}
	541
6c41fc8c	542	if (fMinPtTrackInEmcal > 0 && fTracks && fGeom) {
	543	Bool_t trackInEmcalOk = kFALSE;
	544	Int_t ntracks = fTracks->GetEntries();
	545	for (Int_t i = 0; i < ntracks; i++) {
	546	AliVTrack track = static_cast<AliVTrack>(fTracks->At(i));
	547	if (!AcceptTrack(track))
	548	continue;
	549	if (track->Eta() < fGeom->GetArm1EtaMin() \|\| track->Eta() > fGeom->GetArm1EtaMax() \|\|
	550	track->Phi() < fGeom->GetArm1PhiMin() * TMath::DegToRad() \|\| track->Phi() > fGeom->GetArm1PhiMax() * TMath::DegToRad())
	551	continue;
	552	if (track->Pt() > fMinPtTrackInEmcal) {
	553	trackInEmcalOk = kTRUE;
	554	break;
	555	}
	556	}
	557	if (!trackInEmcalOk)
	558	return kFALSE;
	559	}
	560
6f6dc2ba	561
	562	if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
	563	!(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
	564	!(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
	565	return kFALSE;
	566
	567
2ac8ca4f	568	return kTRUE;
	569	}
	570
c596bd28	571	//________________________________________________________________________
	572	TClonesArray AliAnalysisTaskEmcal::GetArrayFromEvent(const char name, const char *clname)
	573	{
	574	// Get array from event.
	575
	576	TClonesArray *arr = 0;
	577	TString sname(name);
	578	if (!sname.IsNull()) {
	579	arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
	580	if (!arr) {
	581	AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
	582	return 0;
	583	}
af0280a9	584	} else {
af0280a9	585	return 0;
c596bd28	586	}
af0280a9	587
c596bd28	588	if (!clname)
	589	return arr;
	590
	591	TString objname(arr->GetClass()->GetName());
	592	TClass cls(objname);
	593	if (!cls.InheritsFrom(clname)) {
	594	AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
	595	GetName(), cls.GetName(), name, clname));
	596	return 0;
	597	}
	598	return arr;
	599	}
	600
2da09763	601	//________________________________________________________________________
	602	Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
	603	{
	604	// Retrieve objects from event.
	605
2da09763	606	fVertex[0] = 0;
	607	fVertex[1] = 0;
	608	fVertex[2] = 0;
b4339010	609	fNVertCont = 0;
	610
	611	const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
	612	if (vert) {
	613	vert->GetXYZ(fVertex);
	614	fNVertCont = vert->GetNContributors();
	615	}
2da09763	616
1f6fff78	617	fBeamType = GetBeamType();
1f6fff78	618
a52de945	619	if (fBeamType == kAA \|\| fBeamType == kpA ) {
2da09763	620	AliCentrality *aliCent = InputEvent()->GetCentrality();
2da09763	621	if (aliCent) {
493d6e49	622	fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
e304c410	623	if(fNcentBins==4) {
	624	if (fCent >= 0 && fCent < 10) fCentBin = 0;
	625	else if (fCent >= 10 && fCent < 30) fCentBin = 1;
	626	else if (fCent >= 30 && fCent < 50) fCentBin = 2;
	627	else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
	628	else {
	629	AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
	630	fCentBin = 3;
	631	}
	632	}
2da09763	633	else {
e304c410	634	Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
e304c410	635	fCentBin = TMath::FloorNint(fCent/centWidth);
2da09763	636	}
901cc837	637	} else {
901cc837	638	AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
2da09763	639	fCentBin = 3;
2da09763	640	}
3c124212	641	AliEventplane *aliEP = InputEvent()->GetEventplane();
	642	if (aliEP) {
	643	fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
	644	fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
	645	fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
	646	} else {
	647	AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
	648	}
901cc837	649	} else {
2da09763	650	fCent = 99;
	651	fCentBin = 0;
	652	}
	653
2da09763	654	return kTRUE;
2da09763	655	}