[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliEmcalPhysicsSelection.cxx

// $Id$
//
// Emcal physics selection class.
//
// Author: C.Loizides

#include "AliEmcalPhysicsSelection.h"
#include "AliAODEvent.h"
#include "AliESDEvent.h"
#include "AliLog.h"

ClassImp(AliEmcalPhysicsSelection)

AliEmcalPhysicsSelection::AliEmcalPhysicsSelection() : 
  AliPhysicsSelection(), 
  fMarkFastOnly(0), 
  fMarkLedEvent(0),
  fSkipFastOnly(0), 
  fSkipLedEvent(0),
  fCellMinE(-1), 
  fClusMinE(-1),
  fTrackMinPt(-1), 
  fTriggers(0),
  fZvertex(-1),
  fZvertexDiff(0),
  fCentMin(-1),
  fCentMax(-1),
  fIsFastOnly(0),
  fIsLedEvent(0),
  fIsGoodEvent(0),
  fCellMaxE(0), 
  fClusMaxE(0),
  fTrackMaxPt(0)
{
  // Default constructor.
}

//__________________________________________________________________________________________________
UInt_t AliEmcalPhysicsSelection::GetSelectionMask(const TObject* obj) 
{ 
  // Calculate selection mask.
  
  const AliVEvent *ev   = dynamic_cast<const AliVEvent*>(obj);
  if (!ev)
    return 0;

  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();

  UInt_t res = 0;
  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(obj);
  const AliAODEvent *aev = 0;
  if (eev) {
    am->LoadBranch("AliESDHeader.");
    am->LoadBranch("AliESDRun.");
    TString title(eev->GetHeader()->GetTitle());
    if (1&&(title.Length()>0)) {
      res = eev->GetHeader()->GetUniqueID();
      res &= 0x4FFFFFFF;
    } else {
      res = IsCollisionCandidate(eev); 
    }
  } else {
    aev = dynamic_cast<const AliAODEvent*>(obj);
    res = aev->GetHeader()->GetOfflineTrigger();
  }

  // return 0, if 0 found
  if (res==0)
    return 0;

  // return 0, if ptrs are not set
  if ((eev==0) && (aev==0))
    return 0;

  if (fTriggers) { // only process given triggers
    if ((res & fTriggers) == 0)
      return res;
  }

  fIsFastOnly  = kFALSE;
  fIsGoodEvent = kFALSE;
  fIsLedEvent  = kFALSE;
  fCellMaxE    = -1;
  fClusMaxE    = -1;
  fTrackMaxPt  = -1;

  if ((res & AliVEvent::kAnyINT)      || 
      (res & AliVEvent::kSemiCentral) || 
      (res & AliVEvent::kCentral)     || 
      (res & AliVEvent::kEMC1)        || 
      (res & AliVEvent::kEMC7)        || 
      (res & AliVEvent::kEMCEJE)      || 
      (res & AliVEvent::kEMCEGA))
    fIsGoodEvent = kTRUE;
  else {
    return 0;
  }

  if (fZvertexDiff || (fZvertex>0)) {
    Double_t vzPRI = +999;
    Double_t vzSPD = -999;
    const AliVVertex *pv = 0;
    if (eev) {
      am->LoadBranch("PrimaryVertex.");
      pv = eev->GetPrimaryVertex();
    } else {
      pv = aev->GetPrimaryVertex();
    }
    if (pv && pv->GetNContributors()>0) {
      vzPRI = pv->GetZ();
    }
    const AliVVertex *sv = 0;
    if (eev) {
      am->LoadBranch("SPDVertex.");
      sv = eev->GetPrimaryVertexSPD();
    } else {
      sv = aev->GetPrimaryVertexSPD();
    }
    if (sv && sv->GetNContributors()>0) {
      vzSPD = sv->GetZ();
    }
    Double_t  dvertex = TMath::Abs(vzPRI-vzSPD);
    // skip events with dvertex>1mm if requested
    // https://indico.cern.ch/getFile.py/access?contribId=4&resId=0&materialId=slides&confId=189624
    // also check on vertex z if requested
    if (fZvertexDiff && (dvertex>0.1))
      fIsGoodEvent = kFALSE;
    if ((fZvertex>0) && (TMath::Abs(dvertex)>fZvertex))
      fIsGoodEvent = kFALSE;
  }

  if ((fCentMin>-1) && (fCentMax>-1)) {
    Double_t v0mcent = -1;
    AliCentrality *centin = 0;
    if (eev) {
      TTree *tree = am->GetTree();
      if (tree->GetBranch("Centrality."))
        am->LoadBranch("Centrality.");
      centin = dynamic_cast<AliCentrality*>(eev->FindListObject("Centrality"));
    } else {
      centin = const_cast<AliAODEvent*>(aev)->GetCentrality();
    }
    if (centin)
      v0mcent = centin->GetCentralityPercentileUnchecked("V0M");
    if ((v0mcent<fCentMin) || (v0mcent>fCentMax))
      fIsGoodEvent = kFALSE;
  }

  AliVCaloCells *cells   = ev->GetEMCALCells();
  const Short_t nCells   = cells->GetNumberOfCells();
    
  // mark LHC11a fast only partition if requested 
  if (res & AliVEvent::kFastOnly) {
    fIsFastOnly = kTRUE;
    if (fMarkFastOnly||fSkipFastOnly) {
      if (nCells>0) {
        AliFatal(Form("Number of cells %d, even though EMCAL should not be in fast only partition.",nCells));
      }
      fIsGoodEvent = kFALSE;
    }
  }

  if (fCellMinE>0) {
    if (eev)
      am->LoadBranch("EMCALCells.");
    for(Int_t iCell=0; iCell<nCells; ++iCell) {
      Short_t cellId = cells->GetCellNumber(iCell);
      Double_t cellE = cells->GetCellAmplitude(cellId);
      if (cellE>fCellMaxE)
        fCellMaxE = cellE;
    }
  }

  if (fClusMinE>0) {
    if (eev)
      am->LoadBranch("CaloClusters");

    const Int_t nCaloClusters = ev->GetNumberOfCaloClusters();
    for(Int_t iClus = 0; iClus<nCaloClusters; ++iClus) {
      AliVCluster *cl = ev->GetCaloCluster(iClus);
      if (!cl->IsEMCAL()) 
        continue;
      Double_t e = cl->E();
      if (e>fClusMaxE)
        fClusMaxE = e;
    }
  }

  if (fTrackMinPt>0) {
    TClonesArray *trks = 0;
    if (eev) {
      am->LoadBranch("PicoTracks");
      trks = dynamic_cast<TClonesArray*>(eev->FindListObject("PicoTracks"));
      if (!trks) {
        am->LoadBranch("Tracks");
        trks = dynamic_cast<TClonesArray*>(eev->FindListObject("Tracks"));
      }
    } else {
      trks = dynamic_cast<TClonesArray*>(aev->FindListObject("tracks"));
    }
    if (trks) {
      const Int_t Ntracks = trks->GetEntriesFast();
      for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
        AliVTrack *track = static_cast<AliVTrack*>(trks->At(iTracks));
        if (!track)
          continue;
        if (aev) { // a bit ugly since cuts are hard coded for now
          AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
          if (!aodtrack->TestFilterBit(256) && !aodtrack->TestFilterBit(512))
            continue;
        }
        Double_t pt = track->Pt();
        if (pt>fTrackMaxPt)
          fTrackMaxPt = pt;
      }
    }
  }

  // bad cell criterion for LHC11a from 
  // https://indico.cern.ch/materialDisplay.py?contribId=4&materialId=slides&confId=147067
  const Int_t runN = ev->GetRunNumber();
  if ((runN>=144871) && (runN<=146860)) { 

    if (eev)
      am->LoadBranch("EMCALCells.");

    // count cells above threshold
    Int_t nCellCount[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
    for(Int_t iCell=0; iCell<nCells; ++iCell) {
      Short_t cellId = cells->GetCellNumber(iCell);
      Double_t cellE = cells->GetCellAmplitude(cellId);
      Int_t sm       = cellId / (24*48);
      if (cellE>0.1)
        ++nCellCount[sm];
    }

    if (nCellCount[4] > 100)
      fIsLedEvent = kTRUE;
    else {
      if ((runN>=146858) && (runN<=146860)) {
        if ((res&AliVEvent::kMB) && (nCellCount[3]>=21))
          fIsLedEvent = kTRUE;
        else if ((res&AliVEvent::kEMC1) && (nCellCount[3]>=35))
          fIsLedEvent = kTRUE;
      }
    }
    if (fIsLedEvent) {
      fIsGoodEvent = kFALSE;
    }
  }

  if (fIsGoodEvent) {
    if (fCellMaxE>fCellMinE)
      res |= kEmcalHC;
    if ((fClusMaxE>fClusMinE) || (fTrackMaxPt>fTrackMinPt))
      res |= kEmcalHT;
    res |= kEmcalOk;
  }

  if ((fSkipLedEvent && fIsLedEvent) ||
      (fSkipFastOnly && fIsFastOnly))
    res = 0;

  return res;
}
Commit	Line	Data
b0a53615	1	// $Id$
b0a53615	2	//
980821ba	3	// Emcal physics selection class.
b0a53615	4	//
cd231d42	5	// Author: C.Loizides
be94a07d	6
be94a07d	7	#include "AliEmcalPhysicsSelection.h"
bb33f819	8	#include "AliAODEvent.h"
be94a07d	9	#include "AliESDEvent.h"
0abef324	10	#include "AliLog.h"
be94a07d	11
	12	ClassImp(AliEmcalPhysicsSelection)
	13
fb81cdb7	14	AliEmcalPhysicsSelection::AliEmcalPhysicsSelection() :
	15	AliPhysicsSelection(),
	16	fMarkFastOnly(0),
	17	fMarkLedEvent(0),
	18	fSkipFastOnly(0),
	19	fSkipLedEvent(0),
2c091b4a	20	fCellMinE(-1),
	21	fClusMinE(-1),
	22	fTrackMinPt(-1),
	23	fTriggers(0),
1dc3ded9	24	fZvertex(-1),
	25	fZvertexDiff(0),
	26	fCentMin(-1),
	27	fCentMax(-1),
fb81cdb7	28	fIsFastOnly(0),
	29	fIsLedEvent(0),
	30	fIsGoodEvent(0),
	31	fCellMaxE(0),
2c091b4a	32	fClusMaxE(0),
2c091b4a	33	fTrackMaxPt(0)
fb81cdb7	34	{
	35	// Default constructor.
	36	}
	37
be94a07d	38	//__________________________________________________________________________________________________
	39	UInt_t AliEmcalPhysicsSelection::GetSelectionMask(const TObject* obj)
	40	{
bb33f819	41	// Calculate selection mask.
	42
	43	const AliVEvent ev = dynamic_cast<const AliVEvent>(obj);
09422fad	44	if (!ev)
bb33f819	45	return 0;
09422fad	46
09422fad	47	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
bb33f819	48
	49	UInt_t res = 0;
	50	const AliESDEvent eev = dynamic_cast<const AliESDEvent>(obj);
	51	const AliAODEvent *aev = 0;
	52	if (eev) {
09422fad	53	am->LoadBranch("AliESDHeader.");
	54	am->LoadBranch("AliESDRun.");
	55	TString title(eev->GetHeader()->GetTitle());
	56	if (1&&(title.Length()>0)) {
	57	res = eev->GetHeader()->GetUniqueID();
	58	res &= 0x4FFFFFFF;
	59	} else {
	60	res = IsCollisionCandidate(eev);
	61	}
bb33f819	62	} else {
	63	aev = dynamic_cast<const AliAODEvent*>(obj);
	64	res = aev->GetHeader()->GetOfflineTrigger();
	65	}
be94a07d	66
1dc3ded9	67	// return 0, if 0 found
f4fa7802	68	if (res==0)
	69	return 0;
	70
1dc3ded9	71	// return 0, if ptrs are not set
	72	if ((eev==0) && (aev==0))
	73	return 0;
	74
2c091b4a	75	if (fTriggers) { // only process given triggers
1dc3ded9	76	if ((res & fTriggers) == 0)
2c091b4a	77	return res;
	78	}
	79
fb81cdb7	80	fIsFastOnly = kFALSE;
	81	fIsGoodEvent = kFALSE;
	82	fIsLedEvent = kFALSE;
f4fa7802	83	fCellMaxE = -1;
	84	fClusMaxE = -1;
	85	fTrackMaxPt = -1;
fb81cdb7	86
1dc3ded9	87	if ((res & AliVEvent::kAnyINT) \|\|
	88	(res & AliVEvent::kSemiCentral) \|\|
	89	(res & AliVEvent::kCentral) \|\|
	90	(res & AliVEvent::kEMC1) \|\|
	91	(res & AliVEvent::kEMC7) \|\|
	92	(res & AliVEvent::kEMCEJE) \|\|
bb33f819	93	(res & AliVEvent::kEMCEGA))
fb81cdb7	94	fIsGoodEvent = kTRUE;
09422fad	95	else {
	96	return 0;
	97	}
fb81cdb7	98
1dc3ded9	99	if (fZvertexDiff \|\| (fZvertex>0)) {
	100	Double_t vzPRI = +999;
	101	Double_t vzSPD = -999;
	102	const AliVVertex *pv = 0;
09422fad	103	if (eev) {
09422fad	104	am->LoadBranch("PrimaryVertex.");
1dc3ded9	105	pv = eev->GetPrimaryVertex();
09422fad	106	} else {
1dc3ded9	107	pv = aev->GetPrimaryVertex();
09422fad	108	}
1dc3ded9	109	if (pv && pv->GetNContributors()>0) {
	110	vzPRI = pv->GetZ();
	111	}
	112	const AliVVertex *sv = 0;
09422fad	113	if (eev) {
09422fad	114	am->LoadBranch("SPDVertex.");
1dc3ded9	115	sv = eev->GetPrimaryVertexSPD();
09422fad	116	} else {
1dc3ded9	117	sv = aev->GetPrimaryVertexSPD();
09422fad	118	}
1dc3ded9	119	if (sv && sv->GetNContributors()>0) {
	120	vzSPD = sv->GetZ();
	121	}
	122	Double_t dvertex = TMath::Abs(vzPRI-vzSPD);
1b91882c	123	// skip events with dvertex>1mm if requested
1dc3ded9	124	// https://indico.cern.ch/getFile.py/access?contribId=4&resId=0&materialId=slides&confId=189624
1dc3ded9	125	// also check on vertex z if requested
1b91882c	126	if (fZvertexDiff && (dvertex>0.1))
1dc3ded9	127	fIsGoodEvent = kFALSE;
	128	if ((fZvertex>0) && (TMath::Abs(dvertex)>fZvertex))
	129	fIsGoodEvent = kFALSE;
	130	}
	131
	132	if ((fCentMin>-1) && (fCentMax>-1)) {
	133	Double_t v0mcent = -1;
	134	AliCentrality *centin = 0;
	135	if (eev) {
6cf0dbe4	136	TTree *tree = am->GetTree();
	137	if (tree->GetBranch("Centrality."))
	138	am->LoadBranch("Centrality.");
1dc3ded9	139	centin = dynamic_cast<AliCentrality*>(eev->FindListObject("Centrality"));
	140	} else {
	141	centin = const_cast<AliAODEvent*>(aev)->GetCentrality();
	142	}
	143	if (centin)
	144	v0mcent = centin->GetCentralityPercentileUnchecked("V0M");
77a933e7	145	if ((v0mcent<fCentMin) \|\| (v0mcent>fCentMax))
1dc3ded9	146	fIsGoodEvent = kFALSE;
	147	}
	148
2c091b4a	149	AliVCaloCells *cells = ev->GetEMCALCells();
fb81cdb7	150	const Short_t nCells = cells->GetNumberOfCells();
	151
	152	// mark LHC11a fast only partition if requested
	153	if (res & AliVEvent::kFastOnly) {
	154	fIsFastOnly = kTRUE;
	155	if (fMarkFastOnly\|\|fSkipFastOnly) {
	156	if (nCells>0) {
	157	AliFatal(Form("Number of cells %d, even though EMCAL should not be in fast only partition.",nCells));
be94a07d	158	}
fb81cdb7	159	fIsGoodEvent = kFALSE;
be94a07d	160	}
	161	}
	162
2c091b4a	163	if (fCellMinE>0) {
	164	if (eev)
	165	am->LoadBranch("EMCALCells.");
	166	for(Int_t iCell=0; iCell<nCells; ++iCell) {
	167	Short_t cellId = cells->GetCellNumber(iCell);
	168	Double_t cellE = cells->GetCellAmplitude(cellId);
	169	if (cellE>fCellMaxE)
	170	fCellMaxE = cellE;
	171	}
fb81cdb7	172	}
2c091b4a	173
	174	if (fClusMinE>0) {
	175	if (eev)
	176	am->LoadBranch("CaloClusters");
	177
	178	const Int_t nCaloClusters = ev->GetNumberOfCaloClusters();
	179	for(Int_t iClus = 0; iClus<nCaloClusters; ++iClus) {
	180	AliVCluster *cl = ev->GetCaloCluster(iClus);
	181	if (!cl->IsEMCAL())
	182	continue;
	183	Double_t e = cl->E();
	184	if (e>fClusMaxE)
	185	fClusMaxE = e;
	186	}
	187	}
	188
	189	if (fTrackMinPt>0) {
	190	TClonesArray *trks = 0;
	191	if (eev) {
	192	am->LoadBranch("PicoTracks");
	193	trks = dynamic_cast<TClonesArray*>(eev->FindListObject("PicoTracks"));
	194	if (!trks) {
	195	am->LoadBranch("Tracks");
	196	trks = dynamic_cast<TClonesArray*>(eev->FindListObject("Tracks"));
	197	}
	198	} else {
	199	trks = dynamic_cast<TClonesArray*>(aev->FindListObject("tracks"));
	200	}
af0280a9	201	if (trks) {
	202	const Int_t Ntracks = trks->GetEntriesFast();
	203	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
	204	AliVTrack track = static_cast<AliVTrack>(trks->At(iTracks));
	205	if (!track)
f4fa7802	206	continue;
af0280a9	207	if (aev) { // a bit ugly since cuts are hard coded for now
	208	AliAODTrack aodtrack = static_cast<AliAODTrack>(track);
	209	if (!aodtrack->TestFilterBit(256) && !aodtrack->TestFilterBit(512))
	210	continue;
	211	}
	212	Double_t pt = track->Pt();
	213	if (pt>fTrackMaxPt)
	214	fTrackMaxPt = pt;
f4fa7802	215	}
2c091b4a	216	}
fb81cdb7	217	}
	218
	219	// bad cell criterion for LHC11a from
	220	// https://indico.cern.ch/materialDisplay.py?contribId=4&materialId=slides&confId=147067
	221	const Int_t runN = ev->GetRunNumber();
	222	if ((runN>=144871) && (runN<=146860)) {
2c091b4a	223
	224	if (eev)
	225	am->LoadBranch("EMCALCells.");
	226
	227	// count cells above threshold
	228	Int_t nCellCount[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
	229	for(Int_t iCell=0; iCell<nCells; ++iCell) {
	230	Short_t cellId = cells->GetCellNumber(iCell);
	231	Double_t cellE = cells->GetCellAmplitude(cellId);
	232	Int_t sm = cellId / (24*48);
	233	if (cellE>0.1)
	234	++nCellCount[sm];
	235	}
	236
fb81cdb7	237	if (nCellCount[4] > 100)
	238	fIsLedEvent = kTRUE;
	239	else {
	240	if ((runN>=146858) && (runN<=146860)) {
	241	if ((res&AliVEvent::kMB) && (nCellCount[3]>=21))
	242	fIsLedEvent = kTRUE;
	243	else if ((res&AliVEvent::kEMC1) && (nCellCount[3]>=35))
	244	fIsLedEvent = kTRUE;
	245	}
	246	}
	247	if (fIsLedEvent) {
	248	fIsGoodEvent = kFALSE;
	249	}
	250	}
	251
77a933e7	252	if (fIsGoodEvent) {
	253	if (fCellMaxE>fCellMinE)
	254	res \|= kEmcalHC;
	255	if ((fClusMaxE>fClusMinE) \|\| (fTrackMaxPt>fTrackMinPt))
	256	res \|= kEmcalHT;
fb81cdb7	257	res \|= kEmcalOk;
77a933e7	258	}
fb81cdb7	259
	260	if ((fSkipLedEvent && fIsLedEvent) \|\|
	261	(fSkipFastOnly && fIsFastOnly))
	262	res = 0;
	263
be94a07d	264	return res;
be94a07d	265	}