[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"

ClassImp(AliAnalysisTaskEmcal)

//________________________________________________________________________
AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() : 
  AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
  fAnaType(kTPC),
  fForceBeamType(kNA),
  fGeneralHistograms(kFALSE),
  fInitialized(kFALSE),
  fCreateHisto(kTRUE),
  fTracksName(),
  fCaloName(),
  fCaloCellsName(),
  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),
  fNcentBins(4),
  fGeom(0),
  fTracks(0),
  fCaloClusters(0),
  fCaloCells(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(),
  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),
  fNcentBins(4),
  fGeom(0),
  fTracks(0),
  fCaloClusters(0),
  fCaloCells(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->GetLabel() > 0) {
    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->IsEMCAL())
    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 (track->GetLabel() != 0) {
    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()) { 
      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;
    }
  }

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