Rename before further modifications
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliAnalysisTaskCheckV0.cxx
1 /**************************************************************************
2  * Author: Boris Hippolyte.                                               *
3  * Contributors are mentioned in the code where appropriate.              *
4  *                                                                        *
5  * Permission to use, copy, modify and distribute this software and its   *
6  * documentation strictly for non-commercial purposes is hereby granted   *
7  * without fee, provided that the above copyright notice appears in all   *
8  * copies and that both the copyright notice and this permission notice   *
9  * appear in the supporting documentation. The authors make no claims     *
10  * about the suitability of this software for any purpose. It is          *
11  * provided "as is" without express or implied warranty.                  *
12  **************************************************************************/
13
14 //-----------------------------------------------------------------
15 //                 AliAnalysisTaskESDCheckV0 class
16 //            This task is for QAing the V0s from the ESD
17 //              Origin: B.H. Nov2007, hippolyt@in2p3.fr
18 //-----------------------------------------------------------------
19 #include "TChain.h"
20 #include "TTree.h"
21 #include "TList.h"
22 #include "TH1F.h"
23 #include "TCanvas.h"
24
25 #include "AliAnalysisTask.h"
26 #include "AliAnalysisManager.h"
27
28 #include "AliESDEvent.h"
29 #include "AliESDInputHandler.h"
30
31 #include "AliESDv0.h"
32
33 #include "AliAnalysisTaskESDCheckV0.h"
34
35 ClassImp(AliAnalysisTaskESDCheckV0)
36
37 //________________________________________________________________________
38 AliAnalysisTaskESDCheckV0::AliAnalysisTaskESDCheckV0(const char *name) 
39   : AliAnalysisTask(name, ""), fESD(0), fListHist(), 
40     fHistTrackMultiplicity(0), fHistV0Multiplicity(0), fHistV0OnFlyStatus(0),
41     fHistV0MultiplicityOff(0), fHistV0Chi2Off(0),
42     fHistDcaV0DaughtersOff(0), fHistV0CosineOfPointingAngleOff(0),
43     fHistV0RadiusOff(0),fHistDcaV0ToPrimVertexOff(0),
44     fHistDcaPosToPrimVertexOff(0),fHistDcaNegToPrimVertexOff(0),
45     fHistMassK0Off(0),fHistMassLambdaOff(0),fHistMassAntiLambdaOff(0),
46     fHistV0MultiplicityOn(0), fHistV0Chi2On(0),
47     fHistDcaV0DaughtersOn(0), fHistV0CosineOfPointingAngleOn(0),
48     fHistV0RadiusOn(0),fHistDcaV0ToPrimVertexOn(0),
49     fHistDcaPosToPrimVertexOn(0),fHistDcaNegToPrimVertexOn(0),
50     fHistMassK0On(0),fHistMassLambdaOn(0),fHistMassAntiLambdaOn(0)
51 {
52   // Constructor
53
54   // Define input and output slots here
55   // Input slot #0 works with a TChain
56   DefineInput(0, TChain::Class());
57   // Output slot #0 writes into a TH1 container
58   DefineOutput(0, TH1F::Class());
59   // Output slot #1 writes into a TList container
60   DefineOutput(1, TList::Class());
61 }
62
63 //________________________________________________________________________
64 void AliAnalysisTaskESDCheckV0::ConnectInputData(Option_t *) 
65 {
66   // Connect ESD or AOD here
67   // Called once
68
69   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
70   if (!tree) {
71     Printf("ERROR: Could not read chain from input slot 0");
72   } else {
73     // Disable all branches and enable only the needed ones
74     // The next two lines are different when data produced as AliESDEvent is read
75     tree->SetBranchStatus("*", kFALSE);
76     tree->SetBranchStatus("fTracks.*", kTRUE);
77     tree->SetBranchStatus("fV0s.*", kTRUE);
78
79     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
80
81     if (!esdH) {
82       Printf("ERROR: Could not get ESDInputHandler");
83     } else
84       fESD = esdH->GetEvent();
85   }
86 }
87
88 //________________________________________________________________________
89 void AliAnalysisTaskESDCheckV0::CreateOutputObjects()
90 {
91   // Create histograms
92   // Called once
93
94   // Distinguish Track and V0 Multiplicity !
95
96   fListHist = new TList();
97
98   if (!fHistTrackMultiplicity) {
99     fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 250, 0, 250);
100     fListHist->Add(fHistTrackMultiplicity);
101   }
102   if (!fHistV0Multiplicity) {
103     fHistV0Multiplicity = new TH1F("fHistV0Multiplicity", "Multiplicity distribution;Number of V0s;Events", 50, 0, 50);
104     fListHist->Add(fHistV0Multiplicity);
105   }
106   if (!fHistV0OnFlyStatus) {
107     fHistV0OnFlyStatus = new TH1F("fHistV0OnFlyStatus", "V0 On fly status;status;Number of V0s", 3, 0, 3);
108     fListHist->Add(fHistV0OnFlyStatus);
109   }
110
111   // V0 offline distributions
112   if (!fHistV0MultiplicityOff) {
113     fHistV0MultiplicityOff = new TH1F("fHistV0MultiplicityOff", "Multiplicity distribution;Number of V0s;Events", 50, 0, 50);
114     fListHist->Add(fHistV0MultiplicityOff);
115   }
116   if (!fHistV0Chi2Off) {
117     fHistV0Chi2Off = new TH1F("fHistV0Chi2Off", "V0 chi2;chi2;Number of V0s", 33, 0, 33);
118     fListHist->Add(fHistV0Chi2Off);
119   }
120   if (!fHistDcaV0DaughtersOff) {
121     fHistDcaV0DaughtersOff = new TH1F("fHistDcaV0DaughtersOff", "DCA between V0 daughters;DCA (cm);Number of V0s", 300, 0, 3);
122     fListHist->Add(fHistDcaV0DaughtersOff);
123   }
124   if (!fHistV0CosineOfPointingAngleOff) {
125     fHistV0CosineOfPointingAngleOff = new TH1F("fHistV0CosineOfPointingAngleOff", "V0 Cosine of Pointing Angle;Number of V0s", 200, 0, 1);
126     fListHist->Add(fHistV0CosineOfPointingAngleOff);
127   }
128   if (!fHistV0RadiusOff) {
129     fHistV0RadiusOff = new TH1F("fHistV0RadiusOff", "V0 decay radius;Radius (cm);Number of V0s", 33, 0, 33);
130     fListHist->Add(fHistV0RadiusOff);
131   }
132   if (!fHistDcaV0ToPrimVertexOff) {
133     fHistDcaV0ToPrimVertexOff = new TH1F("fHistDcaV0ToPrimVertexOff", "DCA of V0 to Prim. Vertex;DCA (cm);Number of V0s", 300, 0, 3);
134     fListHist->Add(fHistDcaV0ToPrimVertexOff);
135   }
136   if (!fHistDcaPosToPrimVertexOff) {
137     fHistDcaPosToPrimVertexOff = new TH1F("fHistDcaPosToPrimVertexOff", "DCA of V0 neg daughter to Prim. Vertex;DCA (cm);Number of V0s", 300, 0, 3);
138     fListHist->Add(fHistDcaPosToPrimVertexOff);
139   }
140   if (!fHistDcaNegToPrimVertexOff) {
141     fHistDcaNegToPrimVertexOff = new TH1F("fHistDcaNegToPrimVertexOff", "DCA of V0 pos daughter to Prim. Vertex;DCA (cm);Number of V0s", 300, 0, 3);
142     fListHist->Add(fHistDcaNegToPrimVertexOff);
143   }
144
145   if (!fHistMassK0Off) {
146     fHistMassK0Off = new TH1F("fMassK0Off","K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts",100,0.4,0.6);
147     fListHist->Add(fHistMassK0Off);
148   }
149   if (!fHistMassLambdaOff) {
150     fHistMassLambdaOff = new TH1F("fMassLambdaOff","#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts",75,1.05,1.2);
151     fListHist->Add(fHistMassLambdaOff);
152   }
153   if (!fHistMassAntiLambdaOff) {
154     fHistMassAntiLambdaOff = new TH1F("fMassAntiLambdaOff","#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts",75,1.05,1.2);
155     fListHist->Add(fHistMassAntiLambdaOff);
156   }
157
158   // V0 on-the-fly distributions
159   if (!fHistV0MultiplicityOn) {
160     fHistV0MultiplicityOn = new TH1F("fHistV0MultiplicityOn", "Multiplicity distribution;Number of V0s;Events", 50, 0, 50);
161     fListHist->Add(fHistV0MultiplicityOn);
162   }
163   if (!fHistV0Chi2On) {
164     fHistV0Chi2On = new TH1F("fHistV0Chi2On", "V0 chi2;chi2;Number of V0s", 33, 0, 33);
165     fListHist->Add(fHistV0Chi2On);
166   }
167   if (!fHistDcaV0DaughtersOn) {
168     fHistDcaV0DaughtersOn = new TH1F("fHistDcaV0DaughtersOn", "DCA between V0 daughters;DCA (cm);Number of V0s", 300, 0, 3);
169     fListHist->Add(fHistDcaV0DaughtersOn);
170   }
171   if (!fHistV0CosineOfPointingAngleOn) {
172     fHistV0CosineOfPointingAngleOn = new TH1F("fHistV0CosineOfPointingAngleOn", "V0 Cosine of Pointing Angle;Number of V0s", 200, 0, 1);
173     fListHist->Add(fHistV0CosineOfPointingAngleOn);
174   }
175   if (!fHistV0RadiusOn) {
176     fHistV0RadiusOn = new TH1F("fHistV0RadiusOn", "V0 decay radius;Radius (cm);Number of V0s", 33, 0, 33);
177     fListHist->Add(fHistV0RadiusOn);
178   }
179   if (!fHistDcaV0ToPrimVertexOn) {
180     fHistDcaV0ToPrimVertexOn = new TH1F("fHistDcaV0ToPrimVertexOn", "DCA of V0 to Prim. Vertex;DCA (cm);Number of V0s", 300, 0, 3);
181     fListHist->Add(fHistDcaV0ToPrimVertexOn);
182   }
183   if (!fHistDcaPosToPrimVertexOn) {
184     fHistDcaPosToPrimVertexOn = new TH1F("fHistDcaPosToPrimVertexOn", "DCA of V0 neg daughter to Prim. Vertex;DCA (cm);Number of V0s", 300, 0, 3);
185     fListHist->Add(fHistDcaPosToPrimVertexOn);
186   }
187   if (!fHistDcaNegToPrimVertexOn) {
188     fHistDcaNegToPrimVertexOn = new TH1F("fHistDcaNegToPrimVertexOn", "DCA of V0 pos daughter to Prim. Vertex;DCA (cm);Number of V0s", 300, 0, 3);
189     fListHist->Add(fHistDcaNegToPrimVertexOn);
190   }
191
192   if (!fHistMassK0On) {
193     fHistMassK0On = new TH1F("fMassK0On","K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts",100,0.4,0.6);
194     fListHist->Add(fHistMassK0On);
195   }
196   if (!fHistMassLambdaOn) {
197     fHistMassLambdaOn = new TH1F("fMassLambdaOn","#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts",75,1.05,1.2);
198     fListHist->Add(fHistMassLambdaOn);
199   }
200   if (!fHistMassAntiLambdaOn) {
201     fHistMassAntiLambdaOn = new TH1F("fMassAntiLambdaOn","#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts",75,1.05,1.2);
202     fListHist->Add(fHistMassAntiLambdaOn);
203   }
204
205 }
206
207 //________________________________________________________________________
208 void AliAnalysisTaskESDCheckV0::Exec(Option_t *) 
209 {
210   // Main loop
211   // Called for each event
212
213   if (!fESD) {
214     Printf("ERROR: fESD not available");
215     return;
216   }
217
218   //  Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
219   fHistTrackMultiplicity->Fill(fESD->GetNumberOfTracks());
220
221   Int_t nv0s = 0;
222   nv0s = fESD->GetNumberOfV0s();
223
224   Int_t    lOnFlyStatus = 0, nv0sOn = 0, nv0sOff = 0;
225   Double_t lChi2V0 = 0;
226   Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
227   Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
228   Double_t lV0CosineOfPointingAngle = 0;
229   Double_t lV0Radius = 0;
230   Double_t lInvMassK0 = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
231
232   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) 
233     {// This is the begining of the V0 loop
234       AliESDv0 *v0 = fESD->GetV0(iV0);
235       if (!v0) continue;
236
237       Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]); 
238
239       lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
240
241       UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPindex());
242       UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetNindex());
243
244       Double_t lMomPos[3]; v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
245       Double_t lMomNeg[3]; v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
246
247       AliESDtrack *pTrack=fESD->GetTrack(lKeyPos);
248       AliESDtrack *nTrack=fESD->GetTrack(lKeyNeg);
249       if (!pTrack || !nTrack) {
250         Printf("ERROR: Could not retreive one of the daughter track");
251         continue;
252       }
253
254       Float_t tDcaPosToPrimVertex[2];
255       if(pTrack) pTrack->GetImpactParameters(tDcaPosToPrimVertex[0],tDcaPosToPrimVertex[1]);
256       else { tDcaPosToPrimVertex[0]=999.;  tDcaPosToPrimVertex[1]=999.;}
257       lDcaPosToPrimVertex = TMath::Sqrt(tDcaPosToPrimVertex[0]*tDcaPosToPrimVertex[0]+tDcaPosToPrimVertex[1]*tDcaPosToPrimVertex[1]);
258
259       Float_t tDcaNegToPrimVertex[2];
260       if(nTrack) nTrack->GetImpactParameters(tDcaNegToPrimVertex[0],tDcaNegToPrimVertex[1]);
261       else { tDcaNegToPrimVertex[0]=999.;  tDcaNegToPrimVertex[1]=999.;}
262       lDcaNegToPrimVertex = TMath::Sqrt(tDcaNegToPrimVertex[0]*tDcaNegToPrimVertex[0]+tDcaNegToPrimVertex[1]*tDcaNegToPrimVertex[1]);
263
264
265       lOnFlyStatus = v0->GetOnFlyStatus();
266       lChi2V0 = v0->GetChi2V0();
267       lDcaV0Daughters = v0->GetDcaV0Daughters();
268       lDcaV0ToPrimVertex = v0->GetD();
269       lV0CosineOfPointingAngle = v0->GetV0CosineOfPointingAngle();
270
271       // Getting invariant mass infos directly from ESD
272       v0->ChangeMassHypothesis(310);
273       lInvMassK0 = v0->GetEffMass();
274       v0->ChangeMassHypothesis(3122);
275       lInvMassLambda = v0->GetEffMass();
276       v0->ChangeMassHypothesis(-3122);
277       lInvMassAntiLambda = v0->GetEffMass();
278
279       fHistV0OnFlyStatus->Fill(lOnFlyStatus);
280       if(!lOnFlyStatus){
281         nv0sOff++;
282         fHistV0Chi2Off->Fill(lChi2V0);
283         fHistDcaV0ToPrimVertexOff->Fill(lDcaV0ToPrimVertex);
284         fHistDcaV0DaughtersOff->Fill(lDcaV0Daughters);
285         fHistV0CosineOfPointingAngleOff->Fill(lV0CosineOfPointingAngle);
286
287         fHistV0RadiusOff->Fill(lV0Radius);
288         fHistDcaPosToPrimVertexOff->Fill(lDcaPosToPrimVertex);
289         fHistDcaNegToPrimVertexOff->Fill(lDcaNegToPrimVertex);
290
291         // Filling invariant mass histos for all candidates
292         fHistMassK0Off->Fill(lInvMassK0);
293         fHistMassLambdaOff->Fill(lInvMassLambda);
294         fHistMassAntiLambdaOff->Fill(lInvMassAntiLambda);
295       }
296       else {
297         nv0sOn++;
298         fHistV0Chi2On->Fill(lChi2V0);
299         fHistDcaV0ToPrimVertexOn->Fill(lDcaV0ToPrimVertex);
300         fHistDcaV0DaughtersOn->Fill(lDcaV0Daughters);
301         fHistV0CosineOfPointingAngleOn->Fill(lV0CosineOfPointingAngle);
302
303         fHistV0RadiusOn->Fill(lV0Radius);
304         fHistDcaPosToPrimVertexOn->Fill(lDcaPosToPrimVertex);
305         fHistDcaNegToPrimVertexOn->Fill(lDcaNegToPrimVertex);
306
307         // Filling invariant mass histos for all candidates
308         fHistMassK0On->Fill(lInvMassK0);
309         fHistMassLambdaOn->Fill(lInvMassLambda);
310         fHistMassAntiLambdaOn->Fill(lInvMassAntiLambda);
311       }
312     }// This is the end of the V0 loop
313
314   fHistV0Multiplicity->Fill(fESD->GetNumberOfV0s());
315   fHistV0MultiplicityOff->Fill(nv0sOff);
316   fHistV0MultiplicityOn->Fill(nv0sOn);
317
318   // Post output data.
319   PostData(0, fHistTrackMultiplicity);
320   PostData(1, fListHist);
321 }      
322
323 //________________________________________________________________________
324 void AliAnalysisTaskESDCheckV0::Terminate(Option_t *) 
325 {
326   // Draw result to the screen
327   // Called once at the end of the query
328
329   fHistTrackMultiplicity = dynamic_cast<TH1F*> (GetOutputData(0));
330   if (!fHistTrackMultiplicity) {
331     Printf("ERROR: fHistTrackMultiplicity not available");
332     return;
333   }
334   fHistV0Multiplicity = dynamic_cast<TH1F*> (((TList*)GetOutputData(1))->FindObject("fHistV0Multiplicity"));
335   if (!fHistV0Multiplicity) {
336     Printf("ERROR: fHistV0Multiplicity not available");
337     return;
338   }
339   fHistV0MultiplicityOff = dynamic_cast<TH1F*> (((TList*)GetOutputData(1))->FindObject("fHistV0MultiplicityOff"));
340   if (!fHistV0MultiplicityOff) {
341     Printf("ERROR: fHistV0MultiplicityOff not available");
342     return;
343   }
344   fHistV0MultiplicityOn = dynamic_cast<TH1F*> (((TList*)GetOutputData(1))->FindObject("fHistV0MultiplicityOn"));
345   if (!fHistV0MultiplicityOn) {
346     Printf("ERROR: fHistV0MultiplicityOn not available");
347     return;
348   }
349    
350   TCanvas *c2 = new TCanvas("AliAnalysisTaskESDCheckV0","Ptot",10,10,510,510);
351   c2->cd(1)->SetLogy();
352
353   fHistTrackMultiplicity->SetMarkerStyle(22);
354   fHistTrackMultiplicity->DrawCopy("E");
355   fHistV0Multiplicity->SetMarkerStyle(26);
356   fHistV0Multiplicity->DrawCopy("ESAME");
357   fHistV0MultiplicityOff->SetMarkerStyle(24);
358   fHistV0MultiplicityOff->DrawCopy("ESAME");
359   fHistV0MultiplicityOn->SetMarkerStyle(20);
360   fHistV0MultiplicityOn->DrawCopy("ESAME");
361 }