]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALJetTasks/AliAnalysisTaskSOH.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
change default name
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliAnalysisTaskSOH.cxx
CommitLineData
e9fd1d62 1// $Id$\r
2//\r
3// Simulation EMCal task.\r
4//\r
5// Author: Saehanseul Oh\r
6\r
7#include "AliAnalysisTaskSOH.h"\r
8\r
9#include <TH1F.h>\r
10#include <TH2F.h>\r
11\r
12#include "AliAnalysisManager.h"\r
13#include "AliAnalysisTask.h"\r
14#include "AliEMCALRecoUtils.h"\r
15#include "AliEMCALTrack.h"\r
16#include "AliESDCaloCluster.h"\r
17#include "AliESDEvent.h"\r
18#include "AliESDInputHandler.h"\r
19#include "AliESDtrack.h"\r
20#include "AliESDtrackCuts.h"\r
21#include "AliExternalTrackParam.h"\r
22#include "AliLog.h"\r
23#include "AliMCEvent.h"\r
24\r
25ClassImp(AliAnalysisTaskSOH)\r
26\r
27//________________________________________________________________________\r
28AliAnalysisTaskSOH::AliAnalysisTaskSOH() :\r
29 AliAnalysisTaskSE(), \r
30 fESD(0), \r
31 fMC(0), \r
32 fEsdTrackCuts(0x0),\r
33 fHybridTrackCuts1(0x0),\r
34 fHybridTrackCuts2(0x0),\r
35 fTrackIndices(0x0),\r
36 fOutputList(0x0), \r
37 fHEventStat(0), \r
38 fHTrkEffParGenPt(0), \r
39 fHTrkEffDetGenPt(0), \r
40 fHTrkEffDetRecPt(0), \r
41 fHEOverPVsPt(0x0),\r
42 fHEMCalResponsePion(0x0), \r
43 fHEMCalResponseElec(0x0),\r
44 fHEMCalResponseProton(0x0), \r
45 fHEMCalRecPhiEtaClus(0x0),\r
46 fHEMCalRecPhiEtaTrk(0x0), \r
47 fHClsPhiEta(0x0),\r
48 fHEMCalRecdPhidEta(0x0),\r
49 fHEMCalRecdPhidEtaP(0x0),\r
50 fHEMCalRecdPhidEtaM(0x0),\r
51 fHEMCalRecdPhidEta_Truth(0x0),\r
52 fHEMCalRecdPhidEtaP_Truth(0x0),\r
53 fHEMCalRecdPhidEtaM_Truth(0x0),\r
54 fHEMCalRecdPhidEtaposEta(0x0),\r
55 fHEMCalRecdPhidEtanegEta(0x0)\r
56{\r
57 // Constructor\r
58\r
59 // Output slot #1 writes into a TH1 container\r
60 DefineOutput(1, TList::Class());\r
61}\r
62\r
63\r
64//________________________________________________________________________\r
65AliAnalysisTaskSOH::AliAnalysisTaskSOH(const char *name) :\r
66 AliAnalysisTaskSE(name), \r
67 fESD(0), \r
68 fMC(0), \r
69 fEsdTrackCuts(0x0),\r
70 fHybridTrackCuts1(0x0),\r
71 fHybridTrackCuts2(0x0),\r
72 fTrackIndices(0x0),\r
73 fOutputList(0x0), \r
74 fHEventStat(0), \r
75 fHTrkEffParGenPt(0), \r
76 fHTrkEffDetGenPt(0), \r
77 fHTrkEffDetRecPt(0), \r
78 fHEOverPVsPt(0x0),\r
79 fHEMCalResponsePion(0x0), \r
80 fHEMCalResponseElec(0x0),\r
81 fHEMCalResponseProton(0x0), \r
82 fHEMCalRecPhiEtaClus(0x0),\r
83 fHEMCalRecPhiEtaTrk(0x0), \r
84 fHClsPhiEta(0x0),\r
85 fHEMCalRecdPhidEta(0x0),\r
86 fHEMCalRecdPhidEtaP(0x0),\r
87 fHEMCalRecdPhidEtaM(0x0),\r
88 fHEMCalRecdPhidEta_Truth(0x0),\r
89 fHEMCalRecdPhidEtaP_Truth(0x0),\r
90 fHEMCalRecdPhidEtaM_Truth(0x0),\r
91 fHEMCalRecdPhidEtaposEta(0x0),\r
92 fHEMCalRecdPhidEtanegEta(0x0)\r
93{\r
94 // Constructor\r
95\r
96 // Output slot #1 writes into a TH1 container\r
97 DefineOutput(1, TList::Class());\r
98}\r
99\r
100\r
101//________________________________________________________________________\r
102AliAnalysisTaskSOH::~AliAnalysisTaskSOH()\r
103{\r
104 // Destructor.\r
105\r
106 delete fEsdTrackCuts;\r
107 delete fHybridTrackCuts1;\r
108 delete fHybridTrackCuts2;\r
109 delete fTrackIndices;\r
110}\r
111\r
112\r
113//________________________________________________________________________\r
114void AliAnalysisTaskSOH::UserCreateOutputObjects()\r
115{\r
116 // Create histograms, called once.\r
117\r
118 OpenFile(1);\r
119 \r
120 fOutputList = new TList();\r
121 fOutputList->SetOwner(1);\r
122\r
123 fHEventStat = new TH1F("fHEventStat","Event statistics for analysis",1,0,1);\r
124 fHEventStat->GetXaxis()->SetBinLabel(1,"Event");\r
125 fOutputList->Add(fHEventStat);\r
126\r
127 fHTrkEffParGenPt = new TH1F("fHTrkEffParGenPt", "Particle level truth p_{T} distribution of generated primary charged pions;p_{T}^{gen} (GeV/c)",15,0.15,3.1);\r
128 fOutputList->Add(fHTrkEffParGenPt);\r
129\r
130 fHTrkEffDetGenPt = new TH1F("fHTrkEffDetGenPt", "Detector level truth p_{T} distribution of primary charged pions;p_{T}^{gen} (GeV/c)",15,0.15,3.1);\r
131 fOutputList->Add(fHTrkEffDetGenPt);\r
132\r
133 fHTrkEffDetRecPt = new TH1F("fHTrkEffDetRecPt", "Reconstructed track p_{T} distribution of primary charged pions;p_{T}^{rec} (GeV/c)",15,0.15,3.1);\r
134 fOutputList->Add(fHTrkEffDetRecPt);\r
135\r
136 fHEOverPVsPt = new TH2F("fHEOverPVsPt", "E/P vs track p_{T}; p_{T} (GeV/c); E/P", 20 , 0, 4, 40, 0, 3.2);\r
137 fOutputList->Add(fHEOverPVsPt);\r
138\r
139 fHEMCalResponsePion = new TH2F("fHEMCalResponsePion", "Pion E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);\r
140 fOutputList->Add(fHEMCalResponsePion);\r
141\r
142 fHEMCalResponseElec = new TH2F("fHEMCalResponseElec", "Electron E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);\r
143 fOutputList->Add(fHEMCalResponseElec);\r
144\r
145 fHEMCalResponseProton = new TH2F("fHEMCalResponseProton", "Proton E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);\r
146 fOutputList->Add(fHEMCalResponseElec);\r
147\r
148 fHEMCalRecPhiEtaClus = new TH2F("fHEMCalRecPhiEtaClus","EMCAL Cluster#phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);\r
149 fOutputList->Add(fHEMCalRecPhiEtaClus);\r
150\r
151 fHEMCalRecPhiEtaTrk = new TH2F("fHEMCalRecPhiEtaTrk","EMCAL Track #phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);\r
152 fOutputList->Add(fHEMCalRecPhiEtaTrk);\r
153\r
154 fHClsPhiEta = new TH2F("fHClsPhiEta", "cluster #phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);\r
155 fOutputList->Add(fHClsPhiEta);\r
156\r
157 fHEMCalRecdPhidEta = new TH2F("fHEMCalRecdPhidEta","EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
158 fOutputList->Add(fHEMCalRecdPhidEta);\r
159\r
160 fHEMCalRecdPhidEtaP = new TH2F("fHEMCalRecdPhidEtaP","EMCAL Charge+ Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
161 fOutputList->Add(fHEMCalRecdPhidEtaP);\r
162\r
163 fHEMCalRecdPhidEtaM = new TH2F("fHEMCalRecdPhidEtaM","EMCAL Charge- Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
164 fOutputList->Add(fHEMCalRecdPhidEtaM);\r
165\r
166 fHEMCalRecdPhidEta_Truth = new TH2F("fHEMCalRecdPhidEta_Truth","EMCAL Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
167 fOutputList->Add(fHEMCalRecdPhidEta_Truth);\r
168\r
169 fHEMCalRecdPhidEtaP_Truth = new TH2F("fHEMCalRecdPhidEtaP_Truth","EMCAL Charge+ Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
170 fOutputList->Add(fHEMCalRecdPhidEtaP_Truth);\r
171\r
172 fHEMCalRecdPhidEtaM_Truth = new TH2F("fHEMCalRecdPhidEtaM_Truth","EMCAL Charge- Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
173 fOutputList->Add(fHEMCalRecdPhidEtaM_Truth);\r
174\r
175 fHEMCalRecdPhidEtaposEta = new TH2F("fHEMCalRecdPhidEtaposEta","(+eta track) EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
176 fOutputList->Add(fHEMCalRecdPhidEtaposEta);\r
177\r
178 fHEMCalRecdPhidEtanegEta = new TH2F("fHEMCalRecdPhidEtanegEta","(-eta track) EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
179 fOutputList->Add(fHEMCalRecdPhidEtanegEta);\r
180 fTrackIndices = new TArrayI();\r
181\r
182 PostData(1, fOutputList);\r
183}\r
184\r
185//________________________________________________________________________\r
186void AliAnalysisTaskSOH::UserExec(Option_t *) \r
187{\r
188 // Main loop, called for each event.\r
189\r
190 // Post output data.\r
191 fESD = dynamic_cast<AliESDEvent*>(InputEvent());\r
192 if (!fESD) {\r
193 printf("ERROR: fESD not available\n");\r
194 return;\r
195 }\r
196\r
197 fMC = MCEvent();\r
198 if (!fMC) {\r
199 printf("ERROR: fMC not available\n");\r
200 return;\r
201 }\r
202\r
203 fHEventStat->Fill(0.5);\r
204\r
205 if(fTrackIndices) \r
206 fTrackIndices->Reset();\r
207\r
208 ProcessTrack();\r
209 ProcessMc();\r
210 \r
211 PostData(1, fOutputList);\r
212} \r
213 \r
214//________________________________________________________________________\r
215void AliAnalysisTaskSOH::ProcessTrack()\r
216{\r
217 // Process track.\r
218\r
219 fTrackIndices->Set(fESD->GetNumberOfTracks());\r
220 AliDebug(3,Form("%s:%d Selecting tracks",(char*)__FILE__,__LINE__));\r
221\r
222 Int_t isMth = 0;\r
223 Int_t nTracks = 0;\r
224\r
225 Float_t ClsPos[3] = {-999,-999,-999};\r
226 Double_t emcTrkpos[3] = {-999,-999,-999};\r
227\r
228 for(Int_t itr=0; itr<fESD->GetNumberOfCaloClusters(); itr++)\r
229 {\r
230 AliESDCaloCluster *cluster = fESD->GetCaloCluster(itr);\r
231 if(!cluster) continue;\r
232 cluster->GetPosition(ClsPos);\r
233 TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);\r
234 fHClsPhiEta->Fill(VClsPos.Eta(), VClsPos.Phi());\r
235 }\r
236\r
237 for(Int_t itr=0; itr<fESD->GetNumberOfTracks(); itr++)\r
238 {\r
239 AliESDtrack *esdtrack = fESD->GetTrack(itr);\r
240 if(!esdtrack)continue;\r
241 AliESDtrack *newTrack = GetAcceptTrack(esdtrack);\r
242 if(!newTrack) continue;\r
243 \r
244 Double_t clsE = -1;\r
245 Int_t clsIndex = newTrack->GetEMCALcluster();\r
246 if(newTrack->GetEMCALcluster()>-1)\r
247 {\r
248 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsIndex);\r
249 if(IsGoodCluster(cluster))\r
250 {\r
251 isMth=1;\r
252 \r
253 cluster->GetPosition(ClsPos);\r
254 TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);\r
255 \r
256 AliEMCALTrack EMCTrk(*newTrack);\r
257 if(!EMCTrk.PropagateToGlobal(ClsPos[0], ClsPos[1], ClsPos[2], 0.0, 0.0)) {\r
258 continue;\r
259 }\r
260 EMCTrk.GetXYZ(emcTrkpos);\r
261 TVector3 VemcTrkPos(emcTrkpos[0],emcTrkpos[1],emcTrkpos[2]);\r
262 \r
263 fHEMCalRecPhiEtaClus->Fill(VClsPos.Eta(), VClsPos.Phi());\r
264 fHEMCalRecPhiEtaTrk->Fill(VemcTrkPos.Eta(), VemcTrkPos.Phi());\r
265 \r
266 Double_t dPhi = VClsPos.Phi() - VemcTrkPos.Phi();\r
267 if (dPhi < -1*TMath::Pi()) dPhi += (2*TMath::Pi());\r
268 else if (dPhi > TMath::Pi()) dPhi -= (2*TMath::Pi());\r
269 \r
270 Double_t dEta = VClsPos.Eta() - VemcTrkPos.Eta();\r
271 \r
272 fHEMCalRecdPhidEta->Fill(dEta, dPhi); \r
273\r
274 if((newTrack->GetLabel())>-1 && (newTrack->GetLabel()) < fMC->GetNumberOfTracks())\r
275 {\r
276 AliVParticle *vParticle = fMC->GetTrack(newTrack->GetLabel());\r
277 if(IsGoodMcParticle(vParticle, newTrack->GetLabel()))\r
278 {\r
279 fHEMCalRecdPhidEta_Truth->Fill(dEta, dPhi);\r
280 if(vParticle->Charge() > 0) fHEMCalRecdPhidEtaP_Truth->Fill(dEta, dPhi);\r
281 if(vParticle->Charge() < 0) fHEMCalRecdPhidEtaM_Truth->Fill(dEta, dPhi);\r
282 }\r
283 }\r
284\r
285 if(esdtrack->Charge() > 0) {fHEMCalRecdPhidEtaP->Fill(dEta, dPhi);}\r
286 if(esdtrack->Charge() < 0) {fHEMCalRecdPhidEtaM->Fill(dEta, dPhi);}\r
287\r
288 if(VemcTrkPos.Eta() > 0) fHEMCalRecdPhidEtaposEta->Fill(dEta, dPhi);\r
289 if(VemcTrkPos.Eta() < 0) fHEMCalRecdPhidEtanegEta->Fill(dEta, dPhi);\r
290 \r
291 clsE = cluster->E();\r
292 if(newTrack->P()>0) fHEOverPVsPt->Fill(newTrack->Pt(),clsE/newTrack->P());\r
293 }\r
294 \r
295 Int_t ipart = newTrack->GetLabel();\r
296 if(ipart>-1 && ipart<fMC->GetNumberOfTracks())\r
297 {\r
298 AliVParticle* vParticle = fMC->GetTrack(ipart);\r
299 if(isMth && vParticle)\r
300 {\r
301 if(TMath::Abs(vParticle->PdgCode())==211)\r
302 {\r
303 fHEMCalResponsePion->Fill(newTrack->Pt(),clsE/newTrack->P());\r
304 }\r
305 if(TMath::Abs(vParticle->PdgCode())==11)\r
306 {\r
307 fHEMCalResponseElec->Fill(newTrack->Pt(),clsE/newTrack->P());\r
308 }\r
309 if(TMath::Abs(vParticle->PdgCode())==2212)\r
310 {\r
311 fHEMCalResponseProton->Fill(newTrack->Pt(),clsE/newTrack->P());\r
312 }\r
313 }\r
314 }\r
315 }\r
316 fTrackIndices->AddAt(itr,nTracks);\r
317 nTracks++;\r
318 }\r
319 fTrackIndices->Set(nTracks);\r
320}\r
321\r
322//________________________________________________________________________\r
323void AliAnalysisTaskSOH::ProcessMc()\r
324{\r
325 // Process MC.\r
326\r
327 for(Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)\r
328 {\r
329 AliVParticle* vParticle = fMC->GetTrack(ipart);\r
330 if(!IsGoodMcParticle(vParticle, ipart)) continue;\r
331 Int_t pdgCode = vParticle->PdgCode();\r
332 \r
333 //tracking effciency\r
334 if(TMath::Abs(vParticle->Eta())<0.9 && TMath::Abs(pdgCode==211))\r
335 {\r
336 fHTrkEffParGenPt->Fill(vParticle->Pt());\r
337 for(Int_t j=0; j<fTrackIndices->GetSize(); j++)\r
338 {\r
339 AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(j));\r
340 if(esdtrack && esdtrack->GetLabel()==ipart)\r
341 {\r
342 fHTrkEffDetGenPt->Fill(vParticle->Pt());\r
343 fHTrkEffDetRecPt->Fill(esdtrack->Pt());\r
344 break;\r
345 }\r
346 }\r
347 }\r
348 } \r
349}\r
350\r
351//________________________________________________________________________\r
352AliESDtrack *AliAnalysisTaskSOH::GetAcceptTrack(AliESDtrack *esdtrack)\r
353{\r
354 // Get accepted track.\r
355\r
356 static AliESDtrack newTrack;\r
357 if(fEsdTrackCuts->AcceptTrack(esdtrack))\r
358 {\r
359 esdtrack->Copy(newTrack);\r
360 newTrack.SetTRDQuality(0);\r
361 }\r
362 else if(fHybridTrackCuts1->AcceptTrack(esdtrack))\r
363 {\r
364 if(esdtrack->GetConstrainedParam())\r
365 {\r
366 esdtrack->Copy(newTrack);\r
367 const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();\r
368 newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());\r
369 newTrack.SetTRDQuality(1); \r
370 }\r
371 else \r
372 return 0x0;\r
373 }\r
374 else if(fHybridTrackCuts2->AcceptTrack(esdtrack))\r
375 {\r
376 if(esdtrack->GetConstrainedParam())\r
377 {\r
378 esdtrack->Copy(newTrack);\r
379 const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();\r
380 newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());\r
381 newTrack.SetTRDQuality(2); \r
382 }\r
383 else \r
384 return 0x0;\r
385 }\r
386 else\r
387 {\r
388 return 0x0;\r
389 }\r
390\r
391 return &newTrack;\r
392}\r
393\r
394//________________________________________________________________________\r
395Bool_t AliAnalysisTaskSOH::IsGoodMcParticle(AliVParticle* vParticle, Int_t ipart)\r
396{\r
397 // Return true if good MC particle.\r
398\r
399 if(!vParticle) return kFALSE;\r
400 if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;\r
401 if (TMath::Abs(vParticle->Eta())>2) return kFALSE;\r
402 if(vParticle->Pt()<0.15) return kFALSE;\r
403 return kTRUE;\r
404}\r
405\r
406\r
407//________________________________________________________________________\r
408Bool_t AliAnalysisTaskSOH::IsGoodCluster(AliESDCaloCluster *cluster)\r
409{\r
410 // Return true if good cluster.\r
411\r
412 if(!cluster) return kFALSE;\r
413 if (!cluster->IsEMCAL()) return kFALSE;\r
414 return kTRUE;\r
415}\r
416\r
417\r
418//________________________________________________________________________\r
419void AliAnalysisTaskSOH::Terminate(Option_t *) \r
420{\r
421 // Terminate analysis.\r
422}\r