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