In AliMUONClusterInfo:
[u/mrichter/AliRoot.git] / PWGUD / FirstPhysics / AliAnalysisTaskFirstPhysics.cxx
1
2 #include "AliAnalysisTaskFirstPhysics.h"
3
4 #include "TH1D.h"
5 #include "TH2D.h"
6 #include "TList.h"
7 #include "TString.h"
8
9 #include "AliAnalysisTaskSE.h"
10 #include "AliAnalysisManager.h"
11 #include "AliESDtrackCuts.h"
12 #include "AliESDEvent.h"
13 #include "AliESDHeader.h"
14 #include "AliMCEvent.h"
15 #include "AliMultiplicity.h"
16
17 #include "AliGenEventHeader.h"
18 #include "AliGenPythiaEventHeader.h"
19 #include "AliGenDPMjetEventHeader.h"
20 #include "AliStack.h"
21
22 ClassImp(AliAnalysisTaskFirstPhysics)
23
24 //________________________________________________________________________
25 AliAnalysisTaskFirstPhysics::AliAnalysisTaskFirstPhysics(const char *name) :
26   AliAnalysisTaskSE(name),
27   fESD(0),
28   fMCEvent(0),
29   fOutput(0),
30   fbReadMC(0),
31   fMCProcessType(kProcND),
32   fTrigger(0),
33   fCutTrackPtMin(0.15),
34   fCutTrackPtMax(100),
35   fCutEta(0.8),
36   fCutVertexZ(10)
37 {
38   DefineOutput(1, TList::Class());
39   for (Int_t i = 0; i < knTrackCuts; i ++) {
40     fTrackCuts[i] = 0;
41   }
42 }
43
44 //________________________________________________________________________
45 AliAnalysisTaskFirstPhysics::~AliAnalysisTaskFirstPhysics()
46 {
47   // Destructor. Clean-up the output list, but not the histograms that are put inside
48   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
49   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
50     delete fOutput;
51   }
52   for (Int_t i = 0; i < knTrackCuts; i ++) {
53     delete fTrackCuts[i];
54     fTrackCuts[i] = 0;
55   }
56   if (fTrigger) {
57     delete fTrigger;
58   }
59 }
60
61 //________________________________________________________________________
62 void AliAnalysisTaskFirstPhysics::PrepareOutputList()
63 {
64   if (fOutput) {
65     AliError("fOutput already initialised.");
66     return;
67   }
68   fOutput = new TList();
69   fOutput->SetOwner();
70   TH1::SetDefaultSumw2(kTRUE);
71   if (fTrigger) {
72     AliError("fTrigger is already initialised.");
73     return;
74   }
75   fTrigger = new AliTriggerAnalysis;
76 }
77
78 //________________________________________________________________________
79 void AliAnalysisTaskFirstPhysics::PrepareDefaultTrackCuts()
80 {
81   // quality cut on ITS+TPC tracks
82   fTrackCuts[kTrackCutQGlo] = new AliESDtrackCuts();
83   // TPC
84   fTrackCuts[kTrackCutQGlo]->SetMinNClustersTPC(70);
85   fTrackCuts[kTrackCutQGlo]->SetMaxChi2PerClusterTPC(4);
86   fTrackCuts[kTrackCutQGlo]->SetAcceptKinkDaughters(kFALSE);
87   fTrackCuts[kTrackCutQGlo]->SetRequireTPCRefit(kTRUE);
88   // ITS
89   fTrackCuts[kTrackCutQGlo]->SetRequireITSRefit(kTRUE);
90   fTrackCuts[kTrackCutQGlo]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
91   fTrackCuts[kTrackCutQGlo]->SetEtaRange(-fCutEta, fCutEta);
92   fTrackCuts[kTrackCutQGlo]->SetPtRange(fCutTrackPtMin, fCutTrackPtMax);
93   // quality cut on ITS_SA tracks (complementary to ITS+TPC)
94   fTrackCuts[kTrackCutQITS] = new AliESDtrackCuts();
95   fTrackCuts[kTrackCutQITS]->SetRequireITSRefit(kTRUE);
96   fTrackCuts[kTrackCutQITS]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
97   fTrackCuts[kTrackCutQITS]->SetEtaRange(-fCutEta, fCutEta);
98   // primary selection for tracks with SPD hits
99   fTrackCuts[kTrackCutDCAwSPD] = new AliESDtrackCuts();
100   fTrackCuts[kTrackCutDCAwSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
101   fTrackCuts[kTrackCutDCAwSPD]->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
102   fTrackCuts[kTrackCutDCAwSPD]->SetMaxDCAToVertexZ(0.5);
103   fTrackCuts[kTrackCutDCAwSPD]->SetEtaRange(-fCutEta, fCutEta);
104   // primary selection for tracks w/o SPD hits
105   fTrackCuts[kTrackCutDCAwoSPD] = new AliESDtrackCuts();
106   fTrackCuts[kTrackCutDCAwoSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
107   fTrackCuts[kTrackCutDCAwoSPD]->SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
108   fTrackCuts[kTrackCutDCAwoSPD]->SetMaxDCAToVertexZ(0.5);
109   fTrackCuts[kTrackCutDCAwoSPD]->SetEtaRange(-fCutEta, fCutEta);
110
111   // tracks without SPD hits
112   fTrackCuts[kTrackCutNoSPD] = new AliESDtrackCuts();
113   fTrackCuts[kTrackCutNoSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
114   // tracks from TPC SA
115   fTrackCuts[kTrackCutTPConly] = new AliESDtrackCuts();
116   fTrackCuts[kTrackCutTPConly]->SetRequireTPCStandAlone(kTRUE);
117 }
118
119 //________________________________________________________________________
120 void AliAnalysisTaskFirstPhysics::UserCreateOutputObjects()
121 {
122   PrepareOutputList();
123   // Define cuts
124   PrepareDefaultTrackCuts();
125
126   // create more histograms here using UserHisto1d and UserHisto2d
127
128   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
129 }
130
131 //________________________________________________________________________
132 bool AliAnalysisTaskFirstPhysics::GetESDEvent()
133 {
134   fESD = 0;
135   // Get a pointer to the reconstructed event
136   AliVEvent *event = InputEvent();
137   if (!event) {
138     AliError("ERROR: Could not retrieve event.");
139     return false;
140   }
141   // try to access the ESD information
142   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
143   if (!esd) {
144     AliError("InputEvent() is not AliESDEvent.");
145     return false;
146   }
147   fESD = esd;
148   return true;
149 }
150
151 //________________________________________________________________________
152 bool AliAnalysisTaskFirstPhysics::CheckVertex()
153 {
154   // check for good reconstructed vertex
155   if (!(fESD->GetPrimaryVertex()->GetStatus()) ||
156       !(fESD->GetPrimaryVertexSPD()->GetStatus())) {
157     return false;
158   }
159   // if vertex is from spd vertexZ, require more stringent cut
160   if (fESD->GetPrimaryVertex()->IsFromVertexerZ() &&
161       (fESD->GetPrimaryVertex()->GetDispersion() > 0.02 ||
162        fESD->GetPrimaryVertex()->GetZRes() > 0.25 )) {
163     return false;
164   }
165   if (TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > fCutVertexZ) {
166     return false;
167   }
168   return true;
169 }
170
171 //________________________________________________________________________
172 bool AliAnalysisTaskFirstPhysics::CheckVertexMC() {
173   if (!fbReadMC) {
174     return false;
175   }
176   if (TMath::Abs(fMCEvent->GetPrimaryVertex()->GetZ()) > fCutVertexZ) {
177     return false;
178   }
179   return true;
180 }
181
182 //________________________________________________________________________
183 bool AliAnalysisTaskFirstPhysics::PrepareMCInfo()
184 {
185   fbReadMC = true;
186   fMCEvent = MCEvent();
187   if (!fMCEvent) {
188     fbReadMC = false;
189   }
190   if (fbReadMC) {
191     AliDebug(4, TString::Format("MC particles: %d", fMCEvent->GetNumberOfTracks()));
192     AliDebug(4, TString::Format("number of particles in the stack: %d", fMCEvent->Stack()->GetNtrack()));
193     AliGenPythiaEventHeader* pythiaGenHeader = 0;
194     AliGenDPMjetEventHeader* dpmHeader = 0;
195     pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fMCEvent->GenEventHeader());
196     dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(fMCEvent->GenEventHeader());
197     // MC event classification
198     Int_t iProcessType = 0;
199     Int_t iSD1flag = 0, iSD2flag = 0, iDDflag = 0,
200       iCDflag = 0, iELflag = 0, iNDflag = 0;
201     if (pythiaGenHeader) {
202       iProcessType = pythiaGenHeader->ProcessType();
203       iSD1flag = 92;
204       iSD2flag = 93;
205       iDDflag = 94;
206       iELflag = 91;
207       iCDflag = -2;
208       iNDflag = -1;
209       AliDebug(5, TString::Format("Pythia event, %d type", iProcessType));
210     } else if (dpmHeader) {
211       iProcessType = dpmHeader->ProcessType();
212       iSD1flag = 5;
213       iSD2flag = 6;
214       iDDflag = 7;
215       iELflag = 2;
216       iCDflag = 4;
217       iNDflag = 1;
218       AliDebug(5, TString::Format("Phojet event, %d type", iProcessType));
219     } else {
220       AliError("neither pythia nor phojet event");
221     }
222     if (iProcessType == iSD1flag) {
223       fMCProcessType = kProcSD1;
224     } else if (iProcessType == iSD2flag) {
225       fMCProcessType = kProcSD2;
226     } else if (iProcessType == iDDflag) {
227       fMCProcessType = kProcDD;
228     } else if (iProcessType == iELflag) {
229       fMCProcessType = kProcEL;
230     } else if (iProcessType == iCDflag) {
231       fMCProcessType = kProcCD;
232     } else if (iProcessType == iNDflag) {
233       fMCProcessType = kProcND;
234     } else {
235       fMCProcessType = kProcIndef;
236     }
237   }
238   return fbReadMC;
239 }
240
241 //________________________________________________________________________
242 void AliAnalysisTaskFirstPhysics::UserExec(Option_t *)
243 {
244   if (!GetESDEvent()) {
245     AliError("Unable to read the ESD");
246     return;
247   }
248   if (!CheckVertex()) {
249     return;
250   }
251
252   // custom code comes here
253
254   PostData(1, fOutput);
255 }
256
257 //________________________________________________________________________
258 TH1D* AliAnalysisTaskFirstPhysics::UserHisto1d(const char *name, const char *title, const char *xlabel, Int_t nbinsx, Double_t xlow, Double_t xup)
259 {
260   TH1D *hist =  new TH1D(name, title, nbinsx, xlow, xup);
261   hist->GetXaxis()->SetTitle(xlabel);
262   fOutput->Add(hist);
263   return hist;
264 }
265
266 //________________________________________________________________________
267 TH2D* AliAnalysisTaskFirstPhysics::UserHisto2d(const char *name, const char *title, const char *xlabel, Int_t nbinsx, Double_t xlow, Double_t xup, const char *ylabel, Int_t nbinsy, Double_t ylow, Double_t yup)
268 {
269   TH2D *hist =  new TH2D(name, title, nbinsx, xlow, xup, nbinsy, ylow, yup);
270   hist->GetXaxis()->SetTitle(xlabel);
271   hist->GetYaxis()->SetTitle(ylabel);
272   fOutput->Add(hist);
273   return hist;
274 }
275
276 //________________________________________________________________________
277 bool AliAnalysisTaskFirstPhysics::GetHisto1FromOutput(const char *name, TH1D *&h) const
278 {
279   h = dynamic_cast<TH1D*> (fOutput->FindObject(name));
280   if (!h) {
281     AliError(TString::Format("Unable to load histogram from output: %s", name));
282     return false;
283   } else {
284     return true;
285   }
286 }
287
288 //________________________________________________________________________
289 bool AliAnalysisTaskFirstPhysics::GetHisto2FromOutput(const char *name, TH2D *&h) const
290 {
291   h = dynamic_cast<TH2D*> (fOutput->FindObject(name));
292   if (!h) {
293     AliError(TString::Format("Unable to load histogram from output: %s", name));
294     return false;
295   } else {
296     return true;
297   }
298 }
299
300 //________________________________________________________________________
301 void AliAnalysisTaskFirstPhysics::Terminate(Option_t *)
302 {
303   fOutput = dynamic_cast<TList*> (GetOutputData(1));
304   if(!fOutput) {
305     AliError("Could not retrieve TList fOutput.");
306     return;
307   }
308 }