[u/mrichter/AliRoot.git] / PWG / EMCAL / 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),
  fMinCellTrackScale(-1),
  fMaxCellTrackScale(-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 = ((AliVAODHeader*)eev->GetHeader())->GetUniqueID();
      res &= 0x4FFFFFFF;
    } else {
      res = IsCollisionCandidate(eev); 
    }
  } else {
    aev = dynamic_cast<const AliAODEvent*>(obj);
    res = ((AliVAODHeader*)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(vzPRI)>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;
    }
  }

  TClonesArray *trks = 0;
  Int_t Ntracks = 0;

  if (fTrackMinPt>0 || fMinCellTrackScale > 0 || fMaxCellTrackScale > 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)
      Ntracks = trks->GetEntriesFast();
  }
  
  Int_t nAccTracks = 0;

  if (fTrackMinPt>0 || fMinCellTrackScale > 0 || fMaxCellTrackScale > 0) {
    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;
      }
      nAccTracks++;
      Double_t pt = track->Pt();
      if (pt>fTrackMaxPt)
	fTrackMaxPt = pt;
    }
  }

  if (fMinCellTrackScale > 0 || fMaxCellTrackScale > 0) {
    if (nCells < fMinCellTrackScale * nAccTracks || nCells > fMaxCellTrackScale * nAccTracks)
      fIsGoodEvent = kFALSE;
  }

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