macro dir
[u/mrichter/AliRoot.git] / 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   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 = eev->GetHeader()->GetUniqueID();
60       res &= 0x4FFFFFFF;
61     } else {
62       res = IsCollisionCandidate(eev); 
63     }
64   } else {
65     aev = dynamic_cast<const AliAODEvent*>(obj);
66     res = 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 }