increment version number
[u/mrichter/AliRoot.git] / PWG4 / omega3pi / AliAnalysisTaskOmegaPi0PiPi.cxx
CommitLineData
f60ddabc 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Boris Polishchuk (Boris.Polishchuk@cern.ch) *
5 * *
6 * Permission to use, copy, modify and distribute this software and its *
7 * documentation strictly for non-commercial purposes is hereby granted *
8 * without fee, provided that the above copyright notice appears in all *
9 * copies and that both the copyright notice and this permission notice *
10 * appear in the supporting documentation. The authors make no claims *
11 * about the suitability of this software for any purpose. It is *
12 * provided "as is" without express or implied warranty. *
13 **************************************************************************/
14
15////////////////////////////////////////////////////////////////////////////
16//------------------------------------------------------------------------//
17// Class used to do omega(782) -> pi0 pi+ pi- analysis. //
18//------------------------------------------------------------------------//
19////////////////////////////////////////////////////////////////////////////
20
21#include "AliAnalysisTaskOmegaPi0PiPi.h"
22#include "TList.h"
23#include "TH1.h"
24#include "TH2.h"
25#include "AliESDEvent.h"
26#include "TRefArray.h"
27#include "TLorentzVector.h"
28#include "AliESDCaloCluster.h"
29#include "AliESDv0.h"
30#include "AliESDtrack.h"
829e60b8 31#include "AliLog.h"
9302ff77 32#include "TString.h"
33#include "AliPHOSGeoUtils.h"
f60ddabc 34
35ClassImp(AliAnalysisTaskOmegaPi0PiPi)
36
37AliAnalysisTaskOmegaPi0PiPi::AliAnalysisTaskOmegaPi0PiPi() :
38 AliAnalysisTaskSE(),fOutputContainer(0x0),fhM2piSel(0x0),fhDxy(0x0),fhMggSel(0x0),
39 fhImpXY(0x0),fhM3pi0to2(0x0),fhM3pi2to4(0x0),fhM3pi4to6(0x0),fhM3pi6to8(0x0),
9302ff77 40 fhM3pi8to10(0x0),fhM3pi10to12(0x0),fhM3pi12(0x0),fhM3pi0to8(0x0),fModules("123")
f60ddabc 41{
42 //Default constructor.
43}
44
45
46AliAnalysisTaskOmegaPi0PiPi::AliAnalysisTaskOmegaPi0PiPi(const char* name) :
47 AliAnalysisTaskSE(name),fOutputContainer(0x0),fhM2piSel(0x0),fhDxy(0x0),fhMggSel(0x0),
48 fhImpXY(0x0),fhM3pi0to2(0x0),fhM3pi2to4(0x0),fhM3pi4to6(0x0),fhM3pi6to8(0x0),
9302ff77 49 fhM3pi8to10(0x0),fhM3pi10to12(0x0),fhM3pi12(0x0),fhM3pi0to8(0x0),fModules("123")
f60ddabc 50{
51 //Constructor used to create a named task.
52
53 DefineOutput(1,TList::Class());
54}
55
56
57AliAnalysisTaskOmegaPi0PiPi::~AliAnalysisTaskOmegaPi0PiPi()
58{
59 //Destructor
60
61 if(fOutputContainer){
62 fOutputContainer->Delete() ;
63 delete fOutputContainer ;
64 }
65}
66
67
68void AliAnalysisTaskOmegaPi0PiPi::UserCreateOutputObjects()
69{
70 //Allocates histograms and puts it to the output container.
71
72 fOutputContainer = new TList();
73
74 fhM2piSel = new TH1F("hM2pi_sel","V0 inv. mass, Dxy<2cm",100,0.,1.);
75 fOutputContainer->Add(fhM2piSel);
76
77 fhDxy = new TH1F("hDxy","Distance from V0 to primary vertex in XY-plane", 500, 0.,500.);
78 fOutputContainer->Add(fhDxy);
79
80 fhMggSel = new TH1F("hMgg_sel","Inv. mass of two clusters, pT>1GeV",100,0.,0.3);
81 fOutputContainer->Add(fhMggSel);
82
83 fhImpXY = new TH1F("hImpXY","Impact parameters in XY-plane",1000,0.,1.);
84 fOutputContainer->Add(fhImpXY);
85
86
87 //M(3pi) vs pT(3pi), ptrack < 2GeV
88 fhM3pi0to2 = new TH2F("hM3pi_0_2","pTrack: 0-2GeV",100,0.,10., 200,0.4,1.);
89 fOutputContainer->Add(fhM3pi0to2);
90
91 //M(3pi) vs pT(3pi), 2GeV < ptrack < 4GeV
92 fhM3pi2to4 = new TH2F("hM3pi_2_4","pTrack: 2-4GeV",100,0.,10., 200,0.4,1.);
93 fOutputContainer->Add(fhM3pi2to4);
94
95 //M(3pi) vs pT(3pi), 4GeV < ptrack < 6GeV
96 fhM3pi4to6 = new TH2F("hM3pi_4_6","pTrack: 4-6GeV",100,0.,10., 200,0.4,1.);
97 fOutputContainer->Add(fhM3pi4to6);
98
99 //M(3pi) vs pT(3pi), 6GeV < ptrack < 8GeV
100 fhM3pi6to8 = new TH2F("hM3pi_6_8","pTrack: 6-8GeV",100,0.,10., 200,0.4,1.);
101 fOutputContainer->Add(fhM3pi6to8);
102
103 //M(3pi) vs pT(3pi), 8GeV < ptrack < 10GeV
104 fhM3pi8to10 = new TH2F("hM3pi_8_10","pTrack: 8-10GeV",100,0.,10., 200,0.4,1.);
105 fOutputContainer->Add(fhM3pi8to10);
106
107 //M(3pi) vs pT(3pi), 10GeV < ptrack < 12GeV
108 fhM3pi10to12 = new TH2F("hM3pi_10_12","pTrack: 10-12GeV",100,0.,10., 200,0.4,1.);
109 fOutputContainer->Add(fhM3pi10to12);
110
111 //M(3pi) vs pT(3pi), ptrack > 12GeV
112 fhM3pi12 = new TH2F("hM3pi_12","pTrack>12GeV",100,0.,10., 200,0.4,1.);
113 fOutputContainer->Add(fhM3pi12);
114
115 //M(3pi) vs pT(3pi), ptrack < 8GeV
116 fhM3pi0to8 = new TH2F("hM3pi_0_8","pTrack: 0-8GeV",100,0.,10., 200,0.4,1.);
117 fOutputContainer->Add(fhM3pi0to8);
118
119}
120
121
122void AliAnalysisTaskOmegaPi0PiPi::UserExec(Option_t* /* option */)
123{
124 //Main function that does all the job.
125
f60ddabc 126 AliVEvent* event = InputEvent();
127 AliESDEvent* esd = (AliESDEvent*)event ;
128
129 Double_t v[3] ; //vertex ;
130 esd->GetVertex()->GetXYZ(v) ;
131 TVector3 vtx(v);
829e60b8 132
133 AliDebug(2,Form("Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z()));
f60ddabc 134
135 TRefArray * caloClustersArr = new TRefArray();
136 esd->GetPHOSClusters(caloClustersArr);
137
138 const Int_t kNumberOfTracks = esd->GetNumberOfTracks();
139 const Int_t kNumberOfV0s = esd->GetNumberOfV0s();
829e60b8 140
141 AliDebug(2,Form("Tracks: %d. V0s: %d\n",kNumberOfTracks,kNumberOfV0s));
f60ddabc 142
143 Float_t xyImp,zImp,imp1,imp2;
144
145 const Int_t kNumberOfPhosClusters = caloClustersArr->GetEntries() ;
829e60b8 146 AliDebug(2,Form("CaloClusters: %d\n", kNumberOfPhosClusters));
93f8dd31 147 if(kNumberOfPhosClusters<2){
148 delete caloClustersArr;
149 return;
150 }
f60ddabc 151 TLorentzVector pc1; //4-momentum of PHOS cluster 1
152 TLorentzVector pc2; //4-momentum of PHOS cluster 2
153 TLorentzVector pc12;
154
155 TLorentzVector ppos; //4-momentum of positive track
156 TLorentzVector pneg; //4-momentum of negative track
157 TLorentzVector ptr12;
158
159 TLorentzVector pomega;
160 Double_t etrack;
161 Double_t p1,p2; // 3-momentum of tracks
9302ff77 162
163 Int_t relid[4] ;
164 UShort_t *absIds;
f60ddabc 165
9302ff77 166 Int_t minCells=3; //Ncells>2
f60ddabc 167 const Double_t kMpi = 0.1396;
9302ff77 168 AliPHOSGeoUtils* geom = new AliPHOSGeoUtils("IHEP","");
169
f60ddabc 170 for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
171 AliESDCaloCluster *c1 = (AliESDCaloCluster *) caloClustersArr->At(iClu);
9302ff77 172 if(c1->GetNCells()<minCells) continue;
173
174 absIds = c1->GetCellsAbsId();
175 geom->AbsToRelNumbering(absIds[0], relid) ;
176 TString phosMod1;
177 phosMod1 += relid[0];
178
179 if(!fModules.Contains(phosMod1.Data())) continue;
f60ddabc 180 c1->GetMomentum(pc1,v);
181
182 for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
183 AliESDCaloCluster *c2 = (AliESDCaloCluster *) caloClustersArr->At(jClu);
184 if(c2->IsEqual(c1)) continue;
9302ff77 185 if(c2->GetNCells()<minCells) continue;
186
187 absIds = c2->GetCellsAbsId();
188 geom->AbsToRelNumbering(absIds[0], relid) ;
189 TString phosMod2;
190 phosMod2 += relid[0];
191
192 if(!fModules.Contains(phosMod2.Data())) continue;
f60ddabc 193 c2->GetMomentum(pc2,v);
194
195 pc12 = pc1+pc2;
829e60b8 196 AliDebug(2,Form("pc12.M(): %.3f\n",pc12.M()));
f60ddabc 197
9302ff77 198 if(pc12.M()<0.100 || pc12.M()>0.160) continue; //not a pi0 candidate!
199 if(pc12.Pt()<1.5) continue; //pT(pi0) > 1.5GeV
f60ddabc 200
201 fhMggSel->Fill(pc12.M());
202
203 for(Int_t iTrack=0; iTrack<kNumberOfTracks; iTrack++) {
204 AliESDtrack* tr1 = esd->GetTrack(iTrack);
205 p1 = tr1->P();
206 tr1->GetImpactParameters(xyImp,zImp);
207 imp1 = TMath::Abs(xyImp);
208 fhImpXY->Fill(imp1);
209 Short_t charge = tr1->Charge();
210 if(imp1>0.004) continue; // not from the primary vertex!
211
212 etrack = TMath::Sqrt(p1*p1 + kMpi*kMpi);
213 ppos.SetPxPyPzE(tr1->Px(),tr1->Py(),tr1->Pz(),etrack);
214
215 for(Int_t jTrack=iTrack; jTrack<kNumberOfTracks; jTrack++) {
216 AliESDtrack* tr2 = esd->GetTrack(jTrack);
217 p2 = tr2->P();
218 if(tr2->Charge()==charge) continue; //same sign track (WRONG!!)
219 tr2->GetImpactParameters(xyImp,zImp);
220 imp2=TMath::Abs(xyImp);
221 if(imp2>0.004) continue; // not from the primary vertex!
222
223 etrack = TMath::Sqrt(p2*p2 + kMpi*kMpi);
224 pneg.SetPxPyPzE(tr2->Px(),tr2->Py(),tr2->Pz(),etrack);
225
226 ptr12 = ppos + pneg;
227
228// printf("ptr12.M()=%.3f, xyImp1=%f, xyImp2=%f, ch1=%d, ch2=%d,
229// p1=%.3f, p2=%.3f, m1=%.3f, m2=%.3f\n",
230// ptr12.M(),imp1,imp2,tr1->Charge(),tr2->Charge(),
231// tr1->P(),tr2->P(),tr1->M(),tr2->M());
232
233 pomega = pc12 + ptr12;
829e60b8 234 AliDebug(2,Form("pomega.M(): %f\n",pomega.M()));
f60ddabc 235
236 if( p1<2. && p2<2. )
237 fhM3pi0to2->Fill(pomega.Pt(),pomega.M());
238
239 if( (p1>=2. && p1<4.) && (p2>=2. && p2<4.) )
240 fhM3pi2to4->Fill(pomega.Pt(),pomega.M());
241
242 if( (p1>=4. && p1<6.) && (p2>=4. && p2<6.) )
243 fhM3pi4to6->Fill(pomega.Pt(),pomega.M());
244
245 if( (p1>=6. && p1<8.) && (p2>=6. && p2<8.) )
246 fhM3pi6to8->Fill(pomega.Pt(),pomega.M());
247
248 if( (p1>=8. && p1<10.) && (p2>=8. && p2<10.) )
249 fhM3pi8to10->Fill(pomega.Pt(),pomega.M());
250
251 if( (p1>=10. && p1<12.) && (p2>=10. && p2<12.) )
252 fhM3pi10to12->Fill(pomega.Pt(),pomega.M());
253
254 if( (p1>=12.) && (p2>=12.) )
255 fhM3pi12->Fill(pomega.Pt(),pomega.M());
256
257 if( p1<8. && p2<8. )
258 fhM3pi0to8->Fill(pomega.Pt(),pomega.M());
259
260 }
261 }
262
263 for(Int_t iVert=0; iVert<kNumberOfV0s; iVert++) {
264 AliESDv0* v0 = esd->GetV0(iVert);
265 AliESDtrack* ptrack = esd->GetTrack(v0->GetPindex()); //positive track
266 AliESDtrack* ntrack = esd->GetTrack(v0->GetNindex()); //negative track
267
268 etrack = TMath::Sqrt(ptrack->P()*ptrack->P() + kMpi*kMpi);
269 ppos.SetPxPyPzE(ptrack->Px(),ptrack->Py(),ptrack->Pz(),etrack);
270
271 etrack = TMath::Sqrt(ntrack->P()*ntrack->P() + kMpi*kMpi);
272 pneg.SetPxPyPzE(ntrack->Px(),ntrack->Py(),ntrack->Pz(),etrack);
273
274 ptr12 = ppos + pneg;
275
276 Double_t dx = vtx.X() - v0->Xv();
277 Double_t dy = vtx.Y() - v0->Yv();
278 Double_t dxy = TMath::Sqrt(dx*dx + dy*dy);
279 fhDxy->Fill(dxy);
829e60b8 280
281 AliDebug(2,Form("V0: dxy=%.3f\n",dxy));
f60ddabc 282
283 if(dxy<2.) fhM2piSel->Fill(ptr12.M());
284 }
285
286 }
287 }
288
289 delete caloClustersArr;
290 PostData(1,fOutputContainer);
291}