new class for omega measurement via 3 pion combination
[u/mrichter/AliRoot.git] / PWG4 / omega3pi / AliAnalysisTaskOmegaPi0PiPi.cxx
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"
31
32 ClassImp(AliAnalysisTaskOmegaPi0PiPi)
33
34 AliAnalysisTaskOmegaPi0PiPi::AliAnalysisTaskOmegaPi0PiPi() :
35   AliAnalysisTaskSE(),fOutputContainer(0x0),fhM2piSel(0x0),fhDxy(0x0),fhMggSel(0x0),
36   fhImpXY(0x0),fhM3pi0to2(0x0),fhM3pi2to4(0x0),fhM3pi4to6(0x0),fhM3pi6to8(0x0),
37   fhM3pi8to10(0x0),fhM3pi10to12(0x0),fhM3pi12(0x0),fhM3pi0to8(0x0)
38 {
39   //Default constructor.
40 }
41
42
43 AliAnalysisTaskOmegaPi0PiPi::AliAnalysisTaskOmegaPi0PiPi(const char* name) :
44   AliAnalysisTaskSE(name),fOutputContainer(0x0),fhM2piSel(0x0),fhDxy(0x0),fhMggSel(0x0),
45   fhImpXY(0x0),fhM3pi0to2(0x0),fhM3pi2to4(0x0),fhM3pi4to6(0x0),fhM3pi6to8(0x0),
46   fhM3pi8to10(0x0),fhM3pi10to12(0x0),fhM3pi12(0x0),fhM3pi0to8(0x0)
47 {
48   //Constructor used to create a named task.
49
50   DefineOutput(1,TList::Class());
51 }
52
53
54 AliAnalysisTaskOmegaPi0PiPi::~AliAnalysisTaskOmegaPi0PiPi()
55 {
56   //Destructor
57   
58   if(fOutputContainer){
59     fOutputContainer->Delete() ; 
60     delete fOutputContainer ;
61   }
62 }
63
64
65 void AliAnalysisTaskOmegaPi0PiPi::UserCreateOutputObjects()
66 {
67   //Allocates histograms and puts it to the output container.
68   
69   fOutputContainer = new TList();
70   
71   fhM2piSel = new TH1F("hM2pi_sel","V0 inv. mass, Dxy<2cm",100,0.,1.);
72   fOutputContainer->Add(fhM2piSel);
73
74   fhDxy = new TH1F("hDxy","Distance from V0 to primary vertex in XY-plane", 500, 0.,500.);
75   fOutputContainer->Add(fhDxy);
76
77   fhMggSel = new TH1F("hMgg_sel","Inv. mass of two clusters, pT>1GeV",100,0.,0.3);
78   fOutputContainer->Add(fhMggSel);
79
80   fhImpXY = new TH1F("hImpXY","Impact parameters in XY-plane",1000,0.,1.);
81   fOutputContainer->Add(fhImpXY);
82
83
84   //M(3pi) vs pT(3pi), ptrack < 2GeV
85   fhM3pi0to2 = new TH2F("hM3pi_0_2","pTrack: 0-2GeV",100,0.,10., 200,0.4,1.);
86   fOutputContainer->Add(fhM3pi0to2);
87
88   //M(3pi) vs pT(3pi), 2GeV < ptrack < 4GeV
89   fhM3pi2to4 = new TH2F("hM3pi_2_4","pTrack: 2-4GeV",100,0.,10., 200,0.4,1.);
90   fOutputContainer->Add(fhM3pi2to4);
91
92   //M(3pi) vs pT(3pi), 4GeV < ptrack < 6GeV
93   fhM3pi4to6 = new TH2F("hM3pi_4_6","pTrack: 4-6GeV",100,0.,10., 200,0.4,1.);
94   fOutputContainer->Add(fhM3pi4to6);
95
96   //M(3pi) vs pT(3pi), 6GeV < ptrack < 8GeV
97   fhM3pi6to8 = new TH2F("hM3pi_6_8","pTrack: 6-8GeV",100,0.,10., 200,0.4,1.);
98   fOutputContainer->Add(fhM3pi6to8);
99
100   //M(3pi) vs pT(3pi), 8GeV < ptrack < 10GeV
101   fhM3pi8to10 = new TH2F("hM3pi_8_10","pTrack: 8-10GeV",100,0.,10., 200,0.4,1.);
102   fOutputContainer->Add(fhM3pi8to10);
103
104   //M(3pi) vs pT(3pi), 10GeV < ptrack < 12GeV
105   fhM3pi10to12 = new TH2F("hM3pi_10_12","pTrack: 10-12GeV",100,0.,10., 200,0.4,1.);
106   fOutputContainer->Add(fhM3pi10to12);
107
108   //M(3pi) vs pT(3pi), ptrack > 12GeV
109   fhM3pi12 = new TH2F("hM3pi_12","pTrack>12GeV",100,0.,10., 200,0.4,1.);
110   fOutputContainer->Add(fhM3pi12);
111
112   //M(3pi) vs pT(3pi), ptrack < 8GeV
113   fhM3pi0to8 = new TH2F("hM3pi_0_8","pTrack: 0-8GeV",100,0.,10., 200,0.4,1.);
114   fOutputContainer->Add(fhM3pi0to8);
115
116 }
117
118
119 void AliAnalysisTaskOmegaPi0PiPi::UserExec(Option_t* /* option */)
120 {
121   //Main function that does all the job.
122   
123   printf("We are within AliAnalysisTaskOmegaPi0PiPi::UserExec.");
124
125   AliVEvent* event = InputEvent();
126   AliESDEvent* esd = (AliESDEvent*)event ;
127
128   Double_t v[3] ; //vertex ;
129   esd->GetVertex()->GetXYZ(v) ;
130   TVector3 vtx(v);
131   printf("Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
132   
133   TRefArray * caloClustersArr  = new TRefArray();
134   esd->GetPHOSClusters(caloClustersArr);
135
136   const Int_t kNumberOfTracks = esd->GetNumberOfTracks();
137   const Int_t kNumberOfV0s = esd->GetNumberOfV0s();
138   printf("Tracks: %d. V0s: %d\n",kNumberOfTracks,kNumberOfV0s);
139   
140   Float_t xyImp,zImp,imp1,imp2;
141   
142   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
143   printf("CaloClusters: %d\n", kNumberOfPhosClusters);
144   if(kNumberOfPhosClusters<2) return;
145   
146   TLorentzVector pc1; //4-momentum of PHOS cluster 1
147   TLorentzVector pc2; //4-momentum of PHOS cluster 2
148   TLorentzVector pc12;
149
150   TLorentzVector ppos; //4-momentum of positive track
151   TLorentzVector pneg; //4-momentum of negative track
152   TLorentzVector ptr12;
153
154   TLorentzVector pomega;
155   Double_t etrack;
156   Double_t p1,p2; // 3-momentum of tracks
157
158   const Double_t kMpi = 0.1396;
159
160   for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
161     AliESDCaloCluster *c1 = (AliESDCaloCluster *) caloClustersArr->At(iClu);
162     c1->GetMomentum(pc1,v);
163
164     for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
165       AliESDCaloCluster *c2 = (AliESDCaloCluster *) caloClustersArr->At(jClu);
166       if(c2->IsEqual(c1)) continue;
167       c2->GetMomentum(pc2,v);
168
169       pc12 = pc1+pc2;
170       printf("pc12.M(): %.3f\n",pc12.M());
171
172       if(pc12.M()<0.115 || pc12.M()>0.155) continue; //not a pi0 candidate!
173       if(pc12.Pt()<1.) continue; //pT(pi0) > 1GeV
174       
175       fhMggSel->Fill(pc12.M());
176       
177       for(Int_t iTrack=0; iTrack<kNumberOfTracks; iTrack++) {
178         AliESDtrack* tr1 = esd->GetTrack(iTrack);
179         p1 = tr1->P();
180         tr1->GetImpactParameters(xyImp,zImp);
181         imp1 = TMath::Abs(xyImp);
182         fhImpXY->Fill(imp1);
183         Short_t charge = tr1->Charge();
184         if(imp1>0.004) continue; // not from the primary vertex!
185
186         etrack = TMath::Sqrt(p1*p1 + kMpi*kMpi);
187         ppos.SetPxPyPzE(tr1->Px(),tr1->Py(),tr1->Pz(),etrack);
188         
189         for(Int_t jTrack=iTrack; jTrack<kNumberOfTracks; jTrack++) {
190           AliESDtrack* tr2 = esd->GetTrack(jTrack);
191           p2 = tr2->P();
192           if(tr2->Charge()==charge) continue; //same sign track (WRONG!!)
193           tr2->GetImpactParameters(xyImp,zImp);
194           imp2=TMath::Abs(xyImp);
195           if(imp2>0.004) continue; // not from the primary vertex!
196           
197           etrack = TMath::Sqrt(p2*p2 + kMpi*kMpi);
198           pneg.SetPxPyPzE(tr2->Px(),tr2->Py(),tr2->Pz(),etrack);
199
200           ptr12 = ppos + pneg;
201
202 //        printf("ptr12.M()=%.3f, xyImp1=%f, xyImp2=%f, ch1=%d, ch2=%d, 
203 //                  p1=%.3f, p2=%.3f, m1=%.3f, m2=%.3f\n",
204 //               ptr12.M(),imp1,imp2,tr1->Charge(),tr2->Charge(),
205 //               tr1->P(),tr2->P(),tr1->M(),tr2->M());
206           
207           pomega = pc12 + ptr12;
208           printf("pomega.M(): %f\n",pomega.M());
209
210           if( p1<2. && p2<2. )
211             fhM3pi0to2->Fill(pomega.Pt(),pomega.M());
212           
213           if( (p1>=2. && p1<4.) && (p2>=2. && p2<4.) ) 
214             fhM3pi2to4->Fill(pomega.Pt(),pomega.M());     
215           
216           if( (p1>=4. && p1<6.) && (p2>=4. && p2<6.) ) 
217             fhM3pi4to6->Fill(pomega.Pt(),pomega.M());
218           
219           if( (p1>=6. && p1<8.) && (p2>=6. && p2<8.) ) 
220             fhM3pi6to8->Fill(pomega.Pt(),pomega.M());
221           
222           if( (p1>=8. && p1<10.) && (p2>=8. && p2<10.) ) 
223             fhM3pi8to10->Fill(pomega.Pt(),pomega.M());
224           
225           if( (p1>=10. && p1<12.) && (p2>=10. && p2<12.) ) 
226             fhM3pi10to12->Fill(pomega.Pt(),pomega.M());
227           
228           if( (p1>=12.) && (p2>=12.) ) 
229             fhM3pi12->Fill(pomega.Pt(),pomega.M());
230           
231           if( p1<8. && p2<8. )
232             fhM3pi0to8->Fill(pomega.Pt(),pomega.M());
233
234         }
235       }
236       
237       for(Int_t iVert=0; iVert<kNumberOfV0s; iVert++) {
238         AliESDv0* v0 = esd->GetV0(iVert);
239         AliESDtrack* ptrack = esd->GetTrack(v0->GetPindex()); //positive track
240         AliESDtrack* ntrack = esd->GetTrack(v0->GetNindex()); //negative track
241
242         etrack = TMath::Sqrt(ptrack->P()*ptrack->P() + kMpi*kMpi);
243         ppos.SetPxPyPzE(ptrack->Px(),ptrack->Py(),ptrack->Pz(),etrack);
244
245         etrack = TMath::Sqrt(ntrack->P()*ntrack->P() + kMpi*kMpi);
246         pneg.SetPxPyPzE(ntrack->Px(),ntrack->Py(),ntrack->Pz(),etrack);
247
248         ptr12 = ppos + pneg;
249
250         Double_t dx = vtx.X() - v0->Xv();
251         Double_t dy = vtx.Y() - v0->Yv();
252         Double_t dxy = TMath::Sqrt(dx*dx + dy*dy);
253         fhDxy->Fill(dxy);
254         printf("V0: dxy=%.3f\n",dxy);
255
256         if(dxy<2.) fhM2piSel->Fill(ptr12.M());
257       } 
258     
259     }
260   }
261   
262   delete caloClustersArr;
263   PostData(1,fOutputContainer);
264 }