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   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
  48
  49   UInt_t res = 0;
  50   const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(obj);
  51   const AliAODEvent *aev = 0;
  52   if (eev) {
  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     }
  62   } else {
  63     aev = dynamic_cast<const AliAODEvent*>(obj);
  64     res = aev->GetHeader()->GetOfflineTrigger();
  65   }
  66
  67   // return 0, if 0 found
  68   if (res==0)
  69     return 0;
  70
  71   // return 0, if ptrs are not set
  72   if ((eev==0) && (aev==0))
  73     return 0;
  74
  75   if (fTriggers) { // only process given triggers
  76     if ((res & fTriggers) == 0)
  77       return res;
  78   }
  79
  80   fIsFastOnly  = kFALSE;
  81   fIsGoodEvent = kFALSE;
  82   fIsLedEvent  = kFALSE;
  83   fCellMaxE    = -1;
  84   fClusMaxE    = -1;
  85   fTrackMaxPt  = -1;
  86
  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)      ||
  93       (res & AliVEvent::kEMCEGA))
  94     fIsGoodEvent = kTRUE;
  95   else {
  96     return 0;
  97   }
  98
  99   if (fZvertexDiff || (fZvertex>0)) {
 100     Double_t vzPRI = +999;
 101     Double_t vzSPD = -999;
 102     const AliVVertex *pv = 0;
 103     if (eev) {
 104       am->LoadBranch("PrimaryVertex.");
 105       pv = eev->GetPrimaryVertex();
 106     } else {
 107       pv = aev->GetPrimaryVertex();
 108     }
 109     if (pv && pv->GetNContributors()>0) {
 110       vzPRI = pv->GetZ();
 111     }
 112     const AliVVertex *sv = 0;
 113     if (eev) {
 114       am->LoadBranch("SPDVertex.");
 115       sv = eev->GetPrimaryVertexSPD();
 116     } else {
 117       sv = aev->GetPrimaryVertexSPD();
 118     }
 119     if (sv && sv->GetNContributors()>0) {
 120       vzSPD = sv->GetZ();
 121     }
 122     Double_t  dvertex = TMath::Abs(vzPRI-vzSPD);
 123     // skip events with dvertex>1mm if requested
 124     // https://indico.cern.ch/getFile.py/access?contribId=4&resId=0&materialId=slides&confId=189624
 125     // also check on vertex z if requested
 126     if (fZvertexDiff && (dvertex>0.1))
 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) {
 136       TTree *tree = am->GetTree();
 137       if (tree->GetBranch("Centrality."))
 138         am->LoadBranch("Centrality.");
 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");
 145     if ((v0mcent<fCentMin) || (v0mcent>fCentMax))
 146       fIsGoodEvent = kFALSE;
 147   }
 148
 149   AliVCaloCells *cells   = ev->GetEMCALCells();
 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));
 158       }
 159       fIsGoodEvent = kFALSE;
 160     }
 161   }
 162
 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     }
 172   }
 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     }
 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)
 206           continue;
 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;
 215       }
 216     }
 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)) {
 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
 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
 252   if (fIsGoodEvent) {
 253     if (fCellMaxE>fCellMinE)
 254       res |= kEmcalHC;
 255     if ((fClusMaxE>fClusMinE) || (fTrackMaxPt>fTrackMinPt))
 256       res |= kEmcalHT;
 257     res |= kEmcalOk;
 258   }
 259
 260   if ((fSkipLedEvent && fIsLedEvent) ||
 261       (fSkipFastOnly && fIsFastOnly))
 262     res = 0;
 263
 264   return res;
 265 }