AliPhysics  9c66e61 (9c66e61)
AliEmcalPhysicsSelection.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Emcal physics selection class.
4 //
5 // Author: C.Loizides
6 
8 #include "AliAODEvent.h"
9 #include "AliESDEvent.h"
10 #include "AliLog.h"
11 
13 
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 //__________________________________________________________________________________________________
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;
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;
270  res |= kEmcalHT;
271  res |= kEmcalOk;
272  }
273 
274  if ((fSkipLedEvent && fIsLedEvent) ||
276  res = 0;
277 
278  return res;
279 }
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
Bool_t fIsGoodEvent
=true if LED event is found
int Int_t
Definition: External.C:63
Double_t fCellMaxE
=true if good EMCAL event
unsigned int UInt_t
Definition: External.C:33
Bool_t fIsLedEvent
=true if FASTONLY event is found
short Short_t
Definition: External.C:23
virtual UInt_t GetSelectionMask(const TObject *obj)
Double_t fClusMaxE
maximum cell energy in event
Double_t fTrackMaxPt
maximum clus energy in event