3 // Simulation EMCal task.
\r
5 // Author: Saehanseul Oh
\r
7 #include "AliAnalysisTaskSOH.h"
\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
23 #include "AliMCEvent.h"
\r
25 ClassImp(AliAnalysisTaskSOH)
\r
27 //________________________________________________________________________
\r
28 AliAnalysisTaskSOH::AliAnalysisTaskSOH() :
\r
29 AliAnalysisTaskSE(),
\r
33 fHybridTrackCuts1(0x0),
\r
34 fHybridTrackCuts2(0x0),
\r
38 fHTrkEffParGenPt(0),
\r
39 fHTrkEffDetGenPt(0),
\r
40 fHTrkEffDetRecPt(0),
\r
42 fHEMCalResponsePion(0x0),
\r
43 fHEMCalResponseElec(0x0),
\r
44 fHEMCalResponseProton(0x0),
\r
45 fHEMCalRecPhiEtaClus(0x0),
\r
46 fHEMCalRecPhiEtaTrk(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
59 // Output slot #1 writes into a TH1 container
\r
60 DefineOutput(1, TList::Class());
\r
64 //________________________________________________________________________
\r
65 AliAnalysisTaskSOH::AliAnalysisTaskSOH(const char *name) :
\r
66 AliAnalysisTaskSE(name),
\r
70 fHybridTrackCuts1(0x0),
\r
71 fHybridTrackCuts2(0x0),
\r
75 fHTrkEffParGenPt(0),
\r
76 fHTrkEffDetGenPt(0),
\r
77 fHTrkEffDetRecPt(0),
\r
79 fHEMCalResponsePion(0x0),
\r
80 fHEMCalResponseElec(0x0),
\r
81 fHEMCalResponseProton(0x0),
\r
82 fHEMCalRecPhiEtaClus(0x0),
\r
83 fHEMCalRecPhiEtaTrk(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
96 // Output slot #1 writes into a TH1 container
\r
97 DefineOutput(1, TList::Class());
\r
101 //________________________________________________________________________
\r
102 AliAnalysisTaskSOH::~AliAnalysisTaskSOH()
\r
106 delete fEsdTrackCuts;
\r
107 delete fHybridTrackCuts1;
\r
108 delete fHybridTrackCuts2;
\r
109 delete fTrackIndices;
\r
113 //________________________________________________________________________
\r
114 void AliAnalysisTaskSOH::UserCreateOutputObjects()
\r
116 // Create histograms, called once.
\r
120 fOutputList = new TList();
\r
121 fOutputList->SetOwner(1);
\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
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
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
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
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
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
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
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
148 fHEMCalRecPhiEtaClus = new TH2F("fHEMCalRecPhiEtaClus","EMCAL Cluster#phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);
\r
149 fOutputList->Add(fHEMCalRecPhiEtaClus);
\r
151 fHEMCalRecPhiEtaTrk = new TH2F("fHEMCalRecPhiEtaTrk","EMCAL Track #phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);
\r
152 fOutputList->Add(fHEMCalRecPhiEtaTrk);
\r
154 fHClsPhiEta = new TH2F("fHClsPhiEta", "cluster #phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);
\r
155 fOutputList->Add(fHClsPhiEta);
\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
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
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
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
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
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
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
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
182 PostData(1, fOutputList);
\r
185 //________________________________________________________________________
\r
186 void AliAnalysisTaskSOH::UserExec(Option_t *)
\r
188 // Main loop, called for each event.
\r
190 // Post output data.
\r
191 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
\r
193 printf("ERROR: fESD not available\n");
\r
199 printf("ERROR: fMC not available\n");
\r
203 fHEventStat->Fill(0.5);
\r
206 fTrackIndices->Reset();
\r
211 PostData(1, fOutputList);
\r
214 //________________________________________________________________________
\r
215 void AliAnalysisTaskSOH::ProcessTrack()
\r
219 fTrackIndices->Set(fESD->GetNumberOfTracks());
\r
220 AliDebug(3,Form("%s:%d Selecting tracks",(char*)__FILE__,__LINE__));
\r
225 Float_t ClsPos[3] = {-999,-999,-999};
\r
226 Double_t emcTrkpos[3] = {-999,-999,-999};
\r
228 for(Int_t itr=0; itr<fESD->GetNumberOfCaloClusters(); itr++)
\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
237 for(Int_t itr=0; itr<fESD->GetNumberOfTracks(); itr++)
\r
239 AliESDtrack *esdtrack = fESD->GetTrack(itr);
\r
240 if(!esdtrack)continue;
\r
241 AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
\r
242 if(!newTrack) continue;
\r
244 Double_t clsE = -1;
\r
245 Int_t clsIndex = newTrack->GetEMCALcluster();
\r
246 if(newTrack->GetEMCALcluster()>-1)
\r
248 AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsIndex);
\r
249 if(IsGoodCluster(cluster))
\r
253 cluster->GetPosition(ClsPos);
\r
254 TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);
\r
256 AliEMCALTrack EMCTrk(*newTrack);
\r
257 if(!EMCTrk.PropagateToGlobal(ClsPos[0], ClsPos[1], ClsPos[2], 0.0, 0.0)) {
\r
260 EMCTrk.GetXYZ(emcTrkpos);
\r
261 TVector3 VemcTrkPos(emcTrkpos[0],emcTrkpos[1],emcTrkpos[2]);
\r
263 fHEMCalRecPhiEtaClus->Fill(VClsPos.Eta(), VClsPos.Phi());
\r
264 fHEMCalRecPhiEtaTrk->Fill(VemcTrkPos.Eta(), VemcTrkPos.Phi());
\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
270 Double_t dEta = VClsPos.Eta() - VemcTrkPos.Eta();
\r
272 fHEMCalRecdPhidEta->Fill(dEta, dPhi);
\r
274 if((newTrack->GetLabel())>-1 && (newTrack->GetLabel()) < fMC->GetNumberOfTracks())
\r
276 AliVParticle *vParticle = fMC->GetTrack(newTrack->GetLabel());
\r
277 if(IsGoodMcParticle(vParticle, newTrack->GetLabel()))
\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
285 if(esdtrack->Charge() > 0) {fHEMCalRecdPhidEtaP->Fill(dEta, dPhi);}
\r
286 if(esdtrack->Charge() < 0) {fHEMCalRecdPhidEtaM->Fill(dEta, dPhi);}
\r
288 if(VemcTrkPos.Eta() > 0) fHEMCalRecdPhidEtaposEta->Fill(dEta, dPhi);
\r
289 if(VemcTrkPos.Eta() < 0) fHEMCalRecdPhidEtanegEta->Fill(dEta, dPhi);
\r
291 clsE = cluster->E();
\r
292 if(newTrack->P()>0) fHEOverPVsPt->Fill(newTrack->Pt(),clsE/newTrack->P());
\r
295 Int_t ipart = newTrack->GetLabel();
\r
296 if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
\r
298 AliVParticle* vParticle = fMC->GetTrack(ipart);
\r
299 if(isMth && vParticle)
\r
301 if(TMath::Abs(vParticle->PdgCode())==211)
\r
303 fHEMCalResponsePion->Fill(newTrack->Pt(),clsE/newTrack->P());
\r
305 if(TMath::Abs(vParticle->PdgCode())==11)
\r
307 fHEMCalResponseElec->Fill(newTrack->Pt(),clsE/newTrack->P());
\r
309 if(TMath::Abs(vParticle->PdgCode())==2212)
\r
311 fHEMCalResponseProton->Fill(newTrack->Pt(),clsE/newTrack->P());
\r
316 fTrackIndices->AddAt(itr,nTracks);
\r
319 fTrackIndices->Set(nTracks);
\r
322 //________________________________________________________________________
\r
323 void AliAnalysisTaskSOH::ProcessMc()
\r
327 for(Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
\r
329 AliVParticle* vParticle = fMC->GetTrack(ipart);
\r
330 if(!IsGoodMcParticle(vParticle, ipart)) continue;
\r
331 Int_t pdgCode = vParticle->PdgCode();
\r
333 //tracking effciency
\r
334 if(TMath::Abs(vParticle->Eta())<0.9 && TMath::Abs(pdgCode==211))
\r
336 fHTrkEffParGenPt->Fill(vParticle->Pt());
\r
337 for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
\r
339 AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(j));
\r
340 if(esdtrack && esdtrack->GetLabel()==ipart)
\r
342 fHTrkEffDetGenPt->Fill(vParticle->Pt());
\r
343 fHTrkEffDetRecPt->Fill(esdtrack->Pt());
\r
351 //________________________________________________________________________
\r
352 AliESDtrack *AliAnalysisTaskSOH::GetAcceptTrack(AliESDtrack *esdtrack)
\r
354 // Get accepted track.
\r
356 static AliESDtrack newTrack;
\r
357 if(fEsdTrackCuts->AcceptTrack(esdtrack))
\r
359 esdtrack->Copy(newTrack);
\r
360 newTrack.SetTRDQuality(0);
\r
362 else if(fHybridTrackCuts1->AcceptTrack(esdtrack))
\r
364 if(esdtrack->GetConstrainedParam())
\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
374 else if(fHybridTrackCuts2->AcceptTrack(esdtrack))
\r
376 if(esdtrack->GetConstrainedParam())
\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
394 //________________________________________________________________________
\r
395 Bool_t AliAnalysisTaskSOH::IsGoodMcParticle(AliVParticle* vParticle, Int_t ipart)
\r
397 // Return true if good MC particle.
\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
407 //________________________________________________________________________
\r
408 Bool_t AliAnalysisTaskSOH::IsGoodCluster(AliESDCaloCluster *cluster)
\r
410 // Return true if good cluster.
\r
412 if(!cluster) return kFALSE;
\r
413 if (!cluster->IsEMCAL()) return kFALSE;
\r
418 //________________________________________________________________________
\r
419 void AliAnalysisTaskSOH::Terminate(Option_t *)
\r
421 // Terminate analysis.
\r