PWG/EMCAL/AliEmcalPhysicsSelection.cxx

   1 // $Id$
   2 //
   3 // Emcal physics selection class.
   4 //
   5 // Author: C.Loizides
   6
   7 #include "AliEmcalPhysicsSelection.h"
   8 #include "AliAODEvent.h"
   9 #include "AliESDEvent.h"
  10 #include "AliLog.h"
  11
  12 ClassImp(AliEmcalPhysicsSelection)
  13
  14 AliEmcalPhysicsSelection::AliEmcalPhysicsSelection() :
  15   AliPhysicsSelection(),
  16   fMarkFastOnly(0),
  17   fMarkLedEvent(0),
  18   fSkipFastOnly(0),
  19   fSkipLedEvent(0),
  20   fCellMinE(-1),
  21   fClusMinE(-1),
  22   fTrackMinPt(-1),
  23   fTriggers(0),
  24   fZvertex(-1),
  25   fZvertexDiff(0),
  26   fCentMin(-1),
  27   fCentMax(-1),
  28   fMinCellTrackScale(-1),
  29   fMaxCellTrackScale(-1),
  30   fIsFastOnly(0),
  31   fIsLedEvent(0),
  32   fIsGoodEvent(0),
  33   fCellMaxE(0),
  34   fClusMaxE(0),
  35   fTrackMaxPt(0)
  36 {
  37   // Default constructor.
  38 }
  39
  40 //__________________________________________________________________________________________________
  41 UInt_t AliEmcalPhysicsSelection::GetSelectionMask(const TObject* obj)
  42 {
  43   // Calculate selection mask.
  44
  45   const AliVEvent *ev   = dynamic_cast<const AliVEvent*>(obj);
  46   if (!ev)
  47     return 0;
  48
  49   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
  50
  51   UInt_t res = 0;
  52   const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(obj);
  53   const AliAODEvent *aev = 0;
  54   if (eev) {
  55     am->LoadBranch("AliESDHeader.");
  56     am->LoadBranch("AliESDRun.");
  57     TString title(eev->GetHeader()->GetTitle());
  58     if (1&&(title.Length()>0)) {
  59       res = ((AliVAODHeader*)eev->GetHeader())->GetUniqueID();
  60       res &= 0x4FFFFFFF;
  61     } else {
  62       res = IsCollisionCandidate(eev);
  63     }
  64   } else {
  65     aev = dynamic_cast<const AliAODEvent*>(obj);
  66     res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
  67   }
  68
  69   // return 0, if 0 found
  70   if (res==0)
  71     return 0;
  72
  73   // return 0, if ptrs are not set
  74   if ((eev==0) && (aev==0))
  75     return 0;
  76
  77   if (fTriggers) { // only process given triggers
  78     if ((res & fTriggers) == 0)
  79       return res;
  80   }
  81
  82   fIsFastOnly  = kFALSE;
  83   fIsGoodEvent = kFALSE;
  84   fIsLedEvent  = kFALSE;
  85   fCellMaxE    = -1;
  86   fClusMaxE    = -1;
  87   fTrackMaxPt  = -1;
  88
  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)      ||
  95       (res & AliVEvent::kEMCEGA))
  96     fIsGoodEvent = kTRUE;
  97   else {
  98     return 0;
  99   }
 100
 101   if (fZvertexDiff || (fZvertex>0)) {
 102     Double_t vzPRI = +999;
 103     Double_t vzSPD = -999;
 104     const AliVVertex *pv = 0;
 105     if (eev) {
 106       am->LoadBranch("PrimaryVertex.");
 107       pv = eev->GetPrimaryVertex();
 108     } else {
 109       pv = aev->GetPrimaryVertex();
 110     }
 111     if (pv && pv->GetNContributors()>0) {
 112       vzPRI = pv->GetZ();
 113     }
 114     const AliVVertex *sv = 0;
 115     if (eev) {
 116       am->LoadBranch("SPDVertex.");
 117       sv = eev->GetPrimaryVertexSPD();
 118     } else {
 119       sv = aev->GetPrimaryVertexSPD();
 120     }
 121     if (sv && sv->GetNContributors()>0) {
 122       vzSPD = sv->GetZ();
 123     }
 124     Double_t  dvertex = TMath::Abs(vzPRI-vzSPD);
 125     // skip events with dvertex>1mm if requested
 126     // https://indico.cern.ch/getFile.py/access?contribId=4&resId=0&materialId=slides&confId=189624
 127     // also check on vertex z if requested
 128     if (fZvertexDiff && (dvertex>0.1))
 129       fIsGoodEvent = kFALSE;
 130     if ((fZvertex>0) && (TMath::Abs(vzPRI)>fZvertex))
 131       fIsGoodEvent = kFALSE;
 132   }
 133
 134   if ((fCentMin>-1) && (fCentMax>-1)) {
 135     Double_t v0mcent = -1;
 136     AliCentrality *centin = 0;
 137     if (eev) {
 138       TTree *tree = am->GetTree();
 139       if (tree->GetBranch("Centrality."))
 140         am->LoadBranch("Centrality.");
 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");
 147     if ((v0mcent<fCentMin) || (v0mcent>fCentMax))
 148       fIsGoodEvent = kFALSE;
 149   }
 150
 151   AliVCaloCells *cells   = ev->GetEMCALCells();
 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));
 160       }
 161       fIsGoodEvent = kFALSE;
 162     }
 163   }
 164
 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     }
 174   }
 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
 191   TClonesArray *trks = 0;
 192   Int_t Ntracks = 0;
 193
 194   if (fTrackMinPt>0 || fMinCellTrackScale > 0 || fMaxCellTrackScale > 0) {
 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     }
 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;
 220       }
 221       nAccTracks++;
 222       Double_t pt = track->Pt();
 223       if (pt>fTrackMaxPt)
 224         fTrackMaxPt = pt;
 225     }
 226   }
 227
 228   if (fMinCellTrackScale > 0 || fMaxCellTrackScale > 0) {
 229     if (nCells < fMinCellTrackScale * nAccTracks || nCells > fMaxCellTrackScale * nAccTracks)
 230       fIsGoodEvent = kFALSE;
 231   }
 232
 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)) {
 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
 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
 266   if (fIsGoodEvent) {
 267     if (fCellMaxE>fCellMinE)
 268       res |= kEmcalHC;
 269     if ((fClusMaxE>fClusMinE) || (fTrackMaxPt>fTrackMinPt))
 270       res |= kEmcalHT;
 271     res |= kEmcalOk;
 272   }
 273
 274   if ((fSkipLedEvent && fIsLedEvent) ||
 275       (fSkipFastOnly && fIsFastOnly))
 276     res = 0;
 277
 278   return res;
 279 }