create general emcal task lib
[u/mrichter/AliRoot.git] / PWG / EMCAL / AliAnalysisTaskEmcal.cxx
1 // $Id: AliAnalysisTaskEmcal.cxx 56756 2012-05-30 05:03:02Z loizides $
2 //
3 // Emcal base analysis task.
4 //
5 // Author: S.Aiola
6
7 #include "AliAnalysisTaskEmcal.h"
8
9 #include <TChain.h>
10 #include <TClonesArray.h>
11 #include <TList.h>
12 #include <TObject.h>
13
14 #include "AliAODEvent.h"
15 #include "AliAnalysisManager.h"
16 #include "AliCentrality.h"
17 #include "AliEMCALGeometry.h"
18 #include "AliESDEvent.h"
19 #include "AliEmcalJet.h"
20 #include "AliEmcalParticle.h"
21 #include "AliEventplane.h"
22 #include "AliInputEventHandler.h"
23 #include "AliLog.h"
24 #include "AliMCParticle.h"
25 #include "AliVCluster.h"
26 #include "AliVEventHandler.h"
27 #include "AliVParticle.h"
28
29 ClassImp(AliAnalysisTaskEmcal)
30
31 //________________________________________________________________________
32 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() : 
33   AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
34   fAnaType(kTPC),
35   fInitialized(kFALSE),
36   fCreateHisto(kTRUE),
37   fTracksName(),
38   fCaloName(),
39   fMinCent(-999),
40   fMaxCent(-999),
41   fMinVz(-999),
42   fMaxVz(-999),
43   fOffTrigger(AliVEvent::kAny),
44   fTrigClass(),
45   fNbins(500),
46   fMinBinPt(0),
47   fMaxBinPt(250),
48   fClusPtCut(0.15),
49   fTrackPtCut(0.15),
50   fMinTrackEta(-0.9),
51   fMaxTrackEta(0.9),
52   fMinTrackPhi(-10),
53   fMaxTrackPhi(10),
54   fClusTimeCutLow(-10),
55   fClusTimeCutUp(10),
56   fTracks(0),
57   fCaloClusters(0),
58   fCent(0),
59   fCentBin(-1),
60   fEPV0(-1.0),
61   fEPV0A(-1.0),
62   fEPV0C(-1.0),
63   fNVertCont(0),
64   fBeamType(kNA),
65   fOutput(0)
66 {
67   // Default constructor.
68
69   fVertex[0] = 0;
70   fVertex[1] = 0;
71   fVertex[2] = 0;
72 }
73
74 //________________________________________________________________________
75 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) : 
76   AliAnalysisTaskSE(name),
77   fAnaType(kTPC),
78   fInitialized(kFALSE),
79   fCreateHisto(histo),
80   fTracksName(),
81   fCaloName(),
82   fMinCent(-999),
83   fMaxCent(-999),
84   fMinVz(-999),
85   fMaxVz(-999),
86   fOffTrigger(AliVEvent::kAny),
87   fTrigClass(),
88   fNbins(500),
89   fMinBinPt(0),
90   fMaxBinPt(250),
91   fClusPtCut(0.15),
92   fTrackPtCut(0.15),
93   fMinTrackEta(-0.9),
94   fMaxTrackEta(0.9),
95   fMinTrackPhi(-10),
96   fMaxTrackPhi(10),
97   fClusTimeCutLow(-10),
98   fClusTimeCutUp(10),
99   fTracks(0),
100   fCaloClusters(0),
101   fCent(0),
102   fCentBin(-1),
103   fEPV0(-1.0),
104   fEPV0A(-1.0),
105   fEPV0C(-1.0),
106   fNVertCont(0),
107   fBeamType(kNA),
108   fOutput(0)
109 {
110   // Standard constructor.
111
112   fVertex[0] = 0;
113   fVertex[1] = 0;
114   fVertex[2] = 0;
115
116   if (fCreateHisto) {
117     DefineOutput(1, TList::Class()); 
118   }
119 }
120
121 //________________________________________________________________________
122 AliAnalysisTaskEmcal::~AliAnalysisTaskEmcal()
123 {
124   // Destructor
125 }
126
127 //________________________________________________________________________
128 void AliAnalysisTaskEmcal::UserExec(Option_t *) 
129 {
130   // Main loop, called for each event.
131
132   if (!fInitialized)
133     ExecOnce();
134
135   if (!fInitialized)
136     return;
137
138   if (!RetrieveEventObjects())
139     return;
140
141   if (!IsEventSelected()) 
142     return;
143
144   if (!Run())
145     return;
146
147   if (!FillHistograms())
148     return;
149     
150   if (fCreateHisto && fOutput) {
151     // information for this iteration of the UserExec in the container
152     PostData(1, fOutput);
153   }
154 }
155
156 //________________________________________________________________________
157 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Bool_t acceptMC) const
158 {
159   // Return true if cluster is accepted.
160
161   if (!clus)
162     return kFALSE;
163
164   if (!clus->IsEMCAL())
165     return kFALSE;
166
167   if (!acceptMC && clus->Chi2() == 100)
168     return kFALSE;
169
170   if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow)
171     return kFALSE;
172
173   TLorentzVector nPart;
174   clus->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
175
176   if (nPart.Et() < fClusPtCut)
177     return kFALSE;
178
179   return kTRUE;
180 }
181
182 //________________________________________________________________________
183 Bool_t AliAnalysisTaskEmcal::AcceptEmcalPart(AliEmcalParticle *part, Bool_t acceptMC) const
184 {
185   // Return true if EMCal particle is accepted.
186
187   if (!part)
188     return kFALSE;
189
190   if (fAnaType == kEMCAL && !part->IsEMCAL())
191     return kFALSE;
192
193   if ((part->IsTrack() && part->Pt() < fTrackPtCut) || (part->IsCluster() && part->Pt() < fClusPtCut))
194     return kFALSE;
195
196   if (!acceptMC && part->IsMC())
197     return kFALSE;
198
199   return kTRUE;
200 }
201
202 //________________________________________________________________________
203 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVTrack *track, Bool_t acceptMC) const
204 {
205   // Return true if track is accepted.
206
207   if (!track)
208     return kFALSE;
209
210   if (!acceptMC && track->GetLabel() == 100)
211     return kFALSE;
212
213   if (track->Pt() < fTrackPtCut)
214     return kFALSE;
215
216   if (track->Eta() < fMinTrackEta || track->Eta() > fMaxTrackEta || 
217       track->Phi() < fMinTrackPhi || track->Phi() > fMaxTrackPhi)
218     return kFALSE;
219   
220   return kTRUE;
221 }
222
223 //________________________________________________________________________
224 void AliAnalysisTaskEmcal::ExecOnce()
225 {
226   // Init the analysis.
227
228   if (!InputEvent()) {
229     AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
230     return;
231   }
232
233   if (!fCaloName.IsNull() && (fAnaType == kEMCAL || fAnaType == kEMCALOnly) && !fCaloClusters) {
234     fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
235     if (!fCaloClusters) {
236       AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCaloName.Data())); 
237       return;
238     } else {
239       TClass *cl = fCaloClusters->GetClass();
240       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
241         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCaloName.Data())); 
242         fCaloClusters = 0;
243         return;
244       }
245     }
246   }
247
248   if (!fTracksName.IsNull() && fAnaType != kEMCALOnly && !fTracks) {
249     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
250     if (!fTracks) {
251       AliError(Form("%s: Could not retrieve tracks %s!", GetName(), fTracksName.Data())); 
252       return;
253     } else {
254       TClass *cl = fTracks->GetClass();
255       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
256         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracksName.Data())); 
257         fTracks = 0;
258         return;
259       }
260     }
261   }
262   SetInitialized();
263 }
264
265 //_____________________________________________________
266 AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
267 {
268   // Get beam type : pp-AA-pA
269   // ESDs have it directly, AODs get it from hardcoded run number ranges
270
271   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
272   if (esd) {
273     const AliESDRun *run = esd->GetESDRun();
274     TString beamType = run->GetBeamType();
275     if (beamType == "p-p")
276       return kpp;
277     else if (beamType == "A-A")
278       return kAA;
279     else if (beamType == "p-A")
280       return kpA;
281     else
282       return kNA;
283   } else {
284     Int_t runNumber = InputEvent()->GetRunNumber();
285     if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
286         (runNumber >= 166529 && runNumber <= 170593))    // LHC11h
287     {
288       return kAA;
289     } else {
290       return kpp;
291     }
292   }  
293 }
294
295 //________________________________________________________________________
296 Bool_t AliAnalysisTaskEmcal::IsEventSelected()
297 {
298   // Check if event is selected
299
300   if (fOffTrigger != AliVEvent::kAny) {
301     UInt_t res = 0;
302     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
303     if (eev) {
304       res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
305     } else {
306       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
307       if (aev) {
308         res = aev->GetHeader()->GetOfflineTrigger();
309       }
310     }
311     if ((res & fOffTrigger) == 0)
312       return kFALSE;
313   }
314
315   if (!fTrigClass.IsNull()) {
316     TString fired;
317     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
318     if (eev) {
319       fired = eev->GetFiredTriggerClasses();
320     } else {
321       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
322       if (aev) {
323         fired = aev->GetFiredTriggerClasses();
324       }
325     }
326     if (!fired.Contains("-B-"))
327       return kFALSE;
328     TObjArray *arr = fTrigClass.Tokenize("|");
329     if (!arr)
330       return kFALSE;
331     Bool_t match = 0;
332     for (Int_t i=0;i<arr->GetEntriesFast();++i) {
333       TObject *obj = arr->At(i);
334       if (!obj)
335         continue;
336       if (fired.Contains(obj->GetName())) {
337         match = 1;
338         break;
339       }
340     }
341     delete arr;
342     if (!match)
343       return kFALSE;
344   }
345
346   if ((fMinCent != -999) && (fMaxCent != -999)) {
347     if (fCent<fMinCent)
348       return kFALSE;
349     if (fCent>fMaxCent)
350       return kFALSE;
351   }
352
353   if ((fMinVz != -999) && (fMaxVz != -999)) {
354     if (fNVertCont == 0 )
355       return kFALSE;
356     Double_t vz = fVertex[2];
357     if (vz<fMinVz)
358       return kFALSE;
359     if (vz>fMaxVz)
360       return kFALSE;
361   }
362
363   return kTRUE;
364 }
365
366 //________________________________________________________________________
367 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
368 {
369   // Get array from event.
370
371   TClonesArray *arr = 0;
372   TString sname(name);
373   if (!sname.IsNull()) {
374     arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
375     if (!arr) {
376       AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name)); 
377       return 0;
378     }
379   } else {
380     return 0;
381   }
382
383   if (!clname)
384     return arr;
385
386   TString objname(arr->GetClass()->GetName());
387   TClass cls(objname);
388   if (!cls.InheritsFrom(clname)) {
389     AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!", 
390                     GetName(), cls.GetName(), name, clname)); 
391     return 0;
392   }
393   return arr;
394 }
395
396 //________________________________________________________________________
397 Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
398 {
399   // Retrieve objects from event.
400
401   fVertex[0] = 0;
402   fVertex[1] = 0;
403   fVertex[2] = 0;
404   fNVertCont = 0;
405
406   const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
407   if (vert) {
408     vert->GetXYZ(fVertex);
409     fNVertCont = vert->GetNContributors();
410   }
411
412   fBeamType = GetBeamType();
413
414   if (fBeamType == kAA) {
415     AliCentrality *aliCent = InputEvent()->GetCentrality();
416     if (aliCent) {
417       fCent = aliCent->GetCentralityPercentile("V0M");
418       if      (fCent >=  0 && fCent <   10) fCentBin = 0;
419       else if (fCent >= 10 && fCent <   30) fCentBin = 1;
420       else if (fCent >= 30 && fCent <   50) fCentBin = 2;
421       else if (fCent >= 50 && fCent <= 100) fCentBin = 3; 
422       else {
423         AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
424         fCentBin = 3;
425       }
426     } else {
427       AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
428       fCentBin = 3;
429     }
430     AliEventplane *aliEP = InputEvent()->GetEventplane();
431     if (aliEP) {
432       fEPV0  = aliEP->GetEventplane("V0" ,InputEvent());
433       fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
434       fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
435     } else {
436       AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
437     }
438   } else {
439     fCent = 99;
440     fCentBin = 0;
441   }
442
443   return kTRUE;
444 }