Update NetParticle: sjena
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFEemcQA.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TH1F.h"
4 #include "TCanvas.h"
5
6 #include "AliAnalysisTask.h"
7 #include "AliAnalysisManager.h"
8
9 #include "AliESDEvent.h"
10 #include "AliESDInputHandler.h"
11 #include "AliESDtrackCuts.h"
12
13 #include "AliAODEvent.h"
14 #include "AliAODHandler.h"
15
16 #include "AliAnalysisTaskHFEemcQA.h"
17
18 //QA task for EMCAL electron analysis 
19
20 ClassImp(AliAnalysisTaskHFEemcQA)
21   //________________________________________________________________________
22   AliAnalysisTaskHFEemcQA::AliAnalysisTaskHFEemcQA(const char *name) 
23 : AliAnalysisTaskSE(name),
24   fVevent(0),
25   fESD(0),
26   fAOD(0),
27   fOutputList(0),
28   fVtxZ(0),
29   fVtxX(0),
30   fVtxY(0),
31   fTrigMulti(0),
32   fHistClustE(0),
33   fEMCClsEtaPhi(0),
34   fNegTrkIDPt(0),
35   fTrkPt(0),
36   fTrketa(0),
37   fTrkphi(0),
38   fdEdx(0),
39   fTPCNpts(0),
40   fHistPtMatch(0),
41   fEMCTrkMatch(0),
42   fEMCTrkPt(0),
43   fEMCTrketa(0),
44   fEMCTrkphi(0),
45   fEMCdEdx(0),
46   fEMCTPCNpts(0),
47   fHistdEdxEop(0),
48   fHistEop(0),
49   fEleCanTPCNpts(0),
50   fEleCanTPCNCls(0),
51   fEleCanITSNCls(0),
52   fEleCanITShit(0),
53   fEleCanSPD1(0),
54   fEleCanSPD2(0),
55   fEleCanSPDBoth(0),
56   fEleCanSPDOr(0)
57 {
58   // Constructor
59
60   // Define input and output slots here
61   // Input slot #0 works with a TChain
62   DefineInput(0, TChain::Class());
63   // Output slot #0 id reserved by the base class for AOD
64   // Output slot #1 writes into a TH1 container
65   DefineOutput(1, TList::Class());
66 }
67
68 //________________________________________________________________________
69 void AliAnalysisTaskHFEemcQA::UserCreateOutputObjects()
70 {
71   // Create histograms
72   // Called once
73
74   fOutputList = new TList();
75   fOutputList->SetOwner();  
76
77   fVtxZ = new TH1F("fVtxZ","Z vertex position;Vtx_{z};counts",1000,-50,50);
78   fOutputList->Add(fVtxZ);
79
80   fVtxY = new TH1F("fVtxY","Y vertex position;Vtx_{y};counts",1000,-50,50);
81   fOutputList->Add(fVtxY);
82
83   fVtxX = new TH1F("fVtxX","X vertex position;Vtx_{x};counts",1000,-50,50);
84   fOutputList->Add(fVtxX);
85
86   fTrigMulti = new TH2F("fTrigMulti","Multiplicity distribution for different triggers; Trigger type; multiplicity",8,-1,7,2000,0,2000);
87   fOutputList->Add(fTrigMulti);
88
89   fHistClustE = new TH1F("fHistClustE", "EMCAL cluster energy distribution; Cluster E;counts", 500, 0.0, 50.0);
90   fOutputList->Add(fHistClustE);
91
92   fEMCClsEtaPhi = new TH2F("fEMCClsEtaPhi","EMCAL cluster #eta and #phi distribution;#eta;#phi",100,-0.9,0.9,200,0,6.3);
93   fOutputList->Add(fEMCClsEtaPhi);
94
95   fNegTrkIDPt = new TH1F("fNegTrkIDPt", "p_{T} distribution of tracks with negative track id;p_{T} (GeV/c);counts", 500, 0.0, 50.0); 
96   fOutputList->Add(fNegTrkIDPt);
97
98   fTrkPt = new TH1F("fTrkPt","p_{T} distribution of all tracks;p_{T} (GeV/c);counts",1000,0,100);
99   fOutputList->Add(fTrkPt);
100
101   fTrketa = new TH1F("fTrketa","All Track #eta distribution;#eta;counts",100,-1.5,1.5);
102   fOutputList->Add(fTrketa);
103
104   fTrkphi = new TH1F("fTrkphi","All Track #phi distribution;#phi;counts",100,0,6.3);
105   fOutputList->Add(fTrkphi);
106
107   fdEdx = new TH2F("fdEdx","All Track dE/dx distribution;p (GeV/c);dE/dx",200,0,20,500,0,100);
108   fOutputList->Add(fdEdx);
109
110   fTPCNpts = new TH2F("fTPCNpts","All track TPC Npoints used for dE/dx calculation;p (GeV/c);N points",200,0,20,200,0.,200.);
111   fOutputList->Add(fTPCNpts);
112
113   fHistPtMatch = new TH1F("fHistPtMatch", "p_{T} distribution of tracks matched to EMCAL;p_{T} (GeV/c);counts",1000, 0.0, 100.0);
114   fOutputList->Add(fHistPtMatch);                                      
115
116   fEMCTrkMatch = new TH2F("fEMCTrkMatch","Distance of EMCAL cluster to its closest track;#phi;z",100,-0.3,0.3,100,-0.3,0.3);
117   fOutputList->Add(fEMCTrkMatch);
118
119   fEMCTrkPt = new TH1F("fEMCTrkPt","p_{T} distribution of tracks with EMCAL cluster;p_{T} (GeV/c);counts",1000,0,100);
120   fOutputList->Add(fEMCTrkPt);
121
122   fEMCTrketa = new TH1F("fEMCTrketa","#eta distribution of tracks matched to EMCAL;#eta;counts",100,-1.5,1.5);
123   fOutputList->Add(fEMCTrketa);
124
125   fEMCTrkphi = new TH1F("fEMCTrkphi","#phi distribution of tracks matched to EMCAL;#phi;counts",100,0,6.3);
126   fOutputList->Add(fEMCTrkphi);
127
128   fEMCdEdx = new TH2F("fEMCdEdx","dE/dx distribution of tracks matched to EMCAL;p (GeV/c);dE/dx",200,0,20,500,0,100);
129   fOutputList->Add(fEMCdEdx);
130
131   fEMCTPCNpts = new TH2F("fEMCTPCNpts","TPC Npoints used for dE/dx for tracks matched to EMCAL;p (GeV/c);N points",200,0,20,200,0.,200.);
132   fOutputList->Add(fEMCTPCNpts);
133
134   fHistEop = new TH2F("fHistEop", "E/p distribution;p_{T} (GeV/c);E/p", 200,0,20,60, 0.0, 3.0);
135   fOutputList->Add(fHistEop);
136
137   fHistdEdxEop = new TH2F("fHistdEdxEop", "E/p vs. dE/dx;E/p;dE/dx", 60, 0.0, 3.0, 500,0,100);
138   fOutputList->Add(fHistdEdxEop);
139
140   fEleCanTPCNpts = new TH2F("fEleCanTPCNpts","TPC Npoints used for dE/dx for electron candidates;p_{T} (GeV/c);N points",200,0,20,200,0,200);
141   fOutputList->Add(fEleCanTPCNpts);
142
143   fEleCanTPCNCls = new TH2F("fEleCanTPCNCls","TPC N clusters for electron candidates;p_{T} (GeV/c);N TPC clusters",200,0,20,171,-0.5,170.5);
144   fOutputList->Add(fEleCanTPCNCls);
145
146   fEleCanITSNCls = new TH2F("fEleCanITSNCls","ITS N clusters for electron candidates;p_{T} (GeV/c);N ITS clusters",200,0,20,8,-0.5,7.5);
147   fOutputList->Add(fEleCanITSNCls);
148
149   fEleCanITShit = new TH1F("fEleCanITShit","ITS hit map;ITS layer;counts",7,-0.5,6.5);
150   fOutputList->Add(fEleCanITShit);
151
152   fEleCanSPD1 = new TH2F("fEleCanSPD1","Hit on SPD layer 1;p_{T} (GeV/c);Hit",200,0,20,1,0,1);
153   fOutputList->Add(fEleCanSPD1);
154
155   fEleCanSPD2 = new TH2F("fEleCanSPD2","Hit on SPD layer 2;p_{T} (GeV/c);Hit",200,0,20,1,0,1);
156   fOutputList->Add(fEleCanSPD2);
157
158   fEleCanSPDBoth = new TH2F("fEleCanSPDBoth","Tracks with hits on both SPD layer;p_{T} (GeV/c);Hit",200,0,20,1,0,1);
159   fOutputList->Add(fEleCanSPDBoth);
160
161   fEleCanSPDOr = new TH2F("fEleCanSPDOr","Tracks with hits on both SPD layer;p_{T} (GeV/c);Hit",200,0,20,1,0,1);
162   fOutputList->Add(fEleCanSPDOr);
163
164   PostData(1,fOutputList);
165 }
166
167 //________________________________________________________________________
168 void AliAnalysisTaskHFEemcQA::UserExec(Option_t *) 
169 {
170   // Main loop
171   // Called for each event
172   // Post output data.
173
174   UInt_t evSelMask=((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
175
176   fVevent = dynamic_cast<AliVEvent*>(InputEvent());
177   if (!fVevent) {
178     printf("ERROR: fVEvent not available\n");
179     return;
180   }
181
182   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
183   if (fESD) {
184     printf("fESD available\n");
185     //return;
186   }
187
188   AliESDtrackCuts* esdTrackCutsH = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
189   esdTrackCutsH->SetMaxDCAToVertexXY(2.4);
190   esdTrackCutsH->SetMaxDCAToVertexZ(3.2);
191   esdTrackCutsH->SetDCAToVertex2D(kTRUE);
192
193   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
194   if (fAOD) {
195     printf("fAOD available\n");
196     //return;
197   }
198
199   int ntracks;
200   ntracks = fVevent->GetNumberOfTracks();
201   printf("There are %d tracks in this event\n",ntracks);
202   //  if(fESD)printf("There are %d tracks in this event\n", fESD->GetNumberOfTracks());
203   //  if(fAOD)printf("There are %d tracks in this event\n", fAOD->GetNumberOfTracks());
204
205   double Zvertex = -100, Xvertex = -100, Yvertex = -100;
206   const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
207   Zvertex = pVtx->GetZ();  
208   Yvertex = pVtx->GetY();  
209   Xvertex = pVtx->GetX();  
210   fVtxZ->Fill(Zvertex);
211   fVtxX->Fill(Xvertex);
212   fVtxY->Fill(Yvertex);
213   /*
214      if(fESD)
215      {
216      const AliESDVertex *pVtx = fESD->GetPrimaryVertex(); 
217      Zvertex = pVtx->GetZ();
218      }
219      if(fAOD)
220      {
221      const AliAODVertex *pVtx = fAOD->GetPrimaryVertex(); 
222      Zvertex = pVtx->GetZ();
223      }
224    */
225   // trigger check
226   
227   TString firedTrigger;
228   TString TriggerEG1("EG1");
229   TString TriggerEG2("EG2");
230   fVevent->GetFiredTriggerClasses();
231   
232   //if(fAOD) firedTrigger = fAOD->GetFiredTriggerClasses();
233   //else if(fESD) firedTrigger = fESD->GetFiredTriggerClasses();
234
235   Bool_t EG1tr = kFALSE;
236   Bool_t EG2tr = kFALSE;
237
238   if(firedTrigger.Contains(TriggerEG1))EG1tr = kTRUE;
239   if(firedTrigger.Contains(TriggerEG2))EG2tr = kTRUE;
240   
241
242   if (fAOD){
243     Double_t multiplicity=fAOD->GetHeader()->GetRefMultiplicity();
244     fTrigMulti->Fill(-0.5, multiplicity);
245     if(evSelMask & AliVEvent::kAny) fTrigMulti->Fill(0.5, multiplicity);
246     if(evSelMask & AliVEvent::kMB) fTrigMulti->Fill(1.5, multiplicity);
247     if(evSelMask & AliVEvent::kINT7) fTrigMulti->Fill(1.5, multiplicity);
248     if(evSelMask & AliVEvent::kINT8) fTrigMulti->Fill(1.5, multiplicity);
249     if(evSelMask & AliVEvent::kEMC1) fTrigMulti->Fill(2.5, multiplicity);
250     if(evSelMask & AliVEvent::kEMC8) fTrigMulti->Fill(3.5, multiplicity);
251     if(evSelMask & AliVEvent::kEMCEJE) fTrigMulti->Fill(4.5, multiplicity);
252     if(evSelMask & AliVEvent::kEMCEGA) fTrigMulti->Fill(5.5, multiplicity);
253     if(evSelMask & AliVEvent::kEMCEGA & EG2tr) fTrigMulti->Fill(5.5, multiplicity);
254     if(evSelMask & AliVEvent::kEMCEGA & EG2tr) fTrigMulti->Fill(5.5, multiplicity);
255   }
256
257   // event selection
258   if(fabs(Zvertex>10.0))return; 
259   
260   //EMCAL cluster information
261   int Nclust = 0;
262   Nclust = fVevent->GetNumberOfCaloClusters();
263   //  if(fESD)Nclust = fESD->GetNumberOfCaloClusters();
264   //  if(fAOD)Nclust = fAOD->GetNumberOfCaloClusters();
265   for(Int_t icl=0; icl<Nclust; icl++)
266   {
267     AliVCluster *clust = 0x0;
268     clust = fVevent->GetCaloCluster(icl);
269     //    if(fESD)clust = (AliVCluster*) fESD->GetCaloCluster(icl);
270     //    if(fAOD)clust = (AliVCluster*) fAOD->GetCaloCluster(icl);
271     if(clust && clust->IsEMCAL())
272     {
273       double clustE = clust->E();
274       float  emcx[3]; // cluster pos
275       clust->GetPosition(emcx);
276       TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
277       double emcphi = clustpos.Phi(); 
278       double emceta = clustpos.Eta();
279       fHistClustE->Fill(clustE);
280       fEMCClsEtaPhi->Fill(emceta,emcphi);
281     }
282   }
283
284   // Track loop 
285   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
286     
287     //---------combine both esd and aod tracks -------
288     AliVParticle* Vtrack = fVevent->GetTrack(iTracks);
289     if (!Vtrack) {
290       printf("ERROR: Could not receive track %d\n", iTracks);
291       continue;
292     }
293     AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack);
294     AliESDtrack *etrack = dynamic_cast<AliESDtrack*>(Vtrack);
295     AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(Vtrack);
296
297     if(fAOD)
298       if(!atrack->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;
299
300     if(fESD)
301       if(!esdTrackCutsH->AcceptTrack(etrack))continue;
302
303     //track properties
304     double dEdx = track->GetTPCsignal();
305     if(track->GetID()<0) fNegTrkIDPt->Fill(track->Pt());
306     fTrkPt->Fill(track->Pt());
307     fTrketa->Fill(track->Eta());
308     fTrkphi->Fill(track->Phi());
309     fdEdx->Fill(track->P(),dEdx);
310     fTPCNpts->Fill(track->P(),track->GetTPCsignalN());
311
312     //track matching to EMCAL
313     int EMCalIndex = -1;
314     EMCalIndex = track->GetEMCALcluster();
315     if(EMCalIndex < 0) continue;
316     fHistPtMatch->Fill(track->Pt());
317
318     AliVCluster *clustMatch = (AliVCluster*)fAOD->GetCaloCluster(EMCalIndex);
319     if(clustMatch && clustMatch->IsEMCAL()) //isnt clustMatch always EMCal cluster?
320     {
321       fEMCTrkMatch->Fill(clustMatch->GetTrackDx(),clustMatch->GetTrackDz());
322       //properties of tracks matched to the EMCAL
323       fEMCTrkPt->Fill(track->Pt());
324       fEMCTrketa->Fill(track->Eta());
325       fEMCTrkphi->Fill(track->Phi());
326       fEMCdEdx->Fill(track->P(),dEdx);
327       fEMCTPCNpts->Fill(track->P(),track->GetTPCsignalN());
328
329       //E/p distribution
330       double clustMatchE = clustMatch->E();
331       double eop = -1.0;
332       if(track->P()>0)eop = clustMatchE/track->P();
333
334       if(track->Pt()>2.0)fHistdEdxEop->Fill(eop,dEdx);
335       fHistEop->Fill(track->Pt(),eop);
336
337       //track properties of EMCAL electron cadidates
338       if(eop>0.8 and eop<1.2){
339         fEleCanTPCNpts->Fill(track->Pt(),track->GetTPCsignalN());
340         fEleCanTPCNCls->Fill(track->Pt(),track->GetTPCNcls());
341
342         Int_t fITSncls=0;
343         for(Int_t l=0;l<6;l++) {
344           if(TESTBIT(track->GetITSClusterMap(),l)) {
345             fEleCanITShit->Fill(l);
346             if(l==0) fEleCanSPD1->Fill(track->Pt(),0.5);
347             if(l==1) fEleCanSPD2->Fill(track->Pt(),0.5);
348             if(l==0 && l==1) fEleCanSPDBoth->Fill(track->Pt(),0.5);
349             if(l==0 || l==1) fEleCanSPDOr->Fill(track->Pt(),0.5);
350             fITSncls++;
351           }
352         }
353         fEleCanITSNCls->Fill(track->Pt(),fITSncls++);
354       }
355     }
356   } //track loop 
357
358   PostData(1, fOutputList);
359 }      
360 //________________________________________________________________________
361 void AliAnalysisTaskHFEemcQA::Terminate(Option_t *) 
362 {
363   // Draw result to the screen
364   // Called once at the end of the query
365
366   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
367   if (!fOutputList) {
368     printf("ERROR: Output list not available\n");
369     return;
370   }
371
372 }