]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaNeutralMeson.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added ptHard to the information in the AliAODMCHeader, filled in the AODHandler
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaNeutralMeson.cxx
CommitLineData
9725fd2a 1/**************************************************************************\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\r
7 * Permission to use, copy, modify and distribute this software and its *\r
8 * documentation strictly for non-commercial purposes is hereby granted *\r
9 * without fee, provided that the above copyright notice appears in all *\r
10 * copies and that both the copyright notice and this permission notice *\r
11 * appear in the supporting documentation. The authors make no claims *\r
12 * about the suitability of this software for any purpose. It is *\r
13 * provided "as is" without express or implied warranty. *\r
14 **************************************************************************/\r
15/* $Id: $ */\r
16\r
17//_________________________________________________________________________\r
18// a general class to fill two-photon and pi0+gamma invariant mass hisograms \r
19// to be used to extract pi0, eta and omega raw yield.\r
20// also for PHOS acceptance\r
21//--\r
22//-- by Renzhuo Wan (Iopp-wuhan) May 13,2009\r
23//_________________________________________________________________________\r
24\r
25// --- ROOT system ---\r
26#include "TH3F.h"\r
27#include "TH2F.h"\r
28//#include "Riostream.h"\r
29#include "TCanvas.h"\r
30#include "TPad.h"\r
31#include "TROOT.h"\r
32#include "TLorentzVector.h"\r
33#include "TString.h"\r
34#include "TParticle.h"\r
35//---- AliRoot system ----\r
36#include "AliAnaNeutralMeson.h"\r
37#include "AliCaloTrackReader.h"\r
38#include "AliCaloPID.h"\r
39#include "AliStack.h"\r
40#include "AliFidutialCut.h"\r
41#include "AliVEvent.h"\r
42#ifdef __PHOSGEO__\r
43 #include "AliPHOSGeoUtils.h"\r
44#endif\r
45\r
46ClassImp(AliAnaNeutralMeson)\r
47\r
48//______________________________________________________________________________\r
49AliAnaNeutralMeson::AliAnaNeutralMeson() : AliAnaPartCorrBaseClass(),\r
50fAnaPi0Eta(0), fAnaOmega(0), fNPID(0),fNmaxMixEv(0), fNAsy(0),fNMod(0), \r
51fNHistos(0), fPerPtBin(0),\r
52fAsyCut(0), fModCut(0),\r
53fDist(0), fAsy(0), fPt(0), fMass(0),\r
54fInvMassCut(0),\r
55fNbinsAsy(0), fMinAsy(0), fMaxAsy(0),\r
56fNbinsPt(0), fMinPt(0), fMaxPt(0),\r
57fNbinsM(0), fMinM(0), fMaxM(0),\r
58fEventsListPi0Eta(0), fEventsListOmega(0),\r
59fPi0Mass(0), fPi0MassPeakWidthCut(0), fhNClusters(0),\r
60fhRecPhoton(0), fhRecPhotonEtaPhi(0),\r
61fReal2Gamma(0), fMix2Gamma(0),fRealOmega(0),\r
62fMixOmegaA(0),fMixOmegaB(0), fMixOmegaC(0),\r
63fRealTwoGammaAsyPtM(0), fMixTwoGammaAsyPtM(0),\r
64fRealPi0GammaAsyPtM(0), fMixAPi0GammaAsyPtM(0), fMixBPi0GammaAsyPtM(0), fMixCPi0GammaAsyPtM(0),\r
65fhPrimPhotonPt(0), fhPrimPhotonAccPt(0), fhPrimPhotonY(0), fhPrimPhotonAccY(0), \r
66fhPrimPhotonPhi(0), fhPrimPhotonAccPhi(0),\r
67fhPrimPi0Pt(0), fhPrimPi0AccPt(0), fhPrimPi0Y(0), fhPrimPi0AccY(0), fhPrimPi0Phi(0), fhPrimPi0AccPhi(0), \r
68fhPrimEtaPt(0), fhPrimEtaAccPt(0), fhPrimEtaY(0), fhPrimEtaAccY(0), fhPrimEtaPhi(0), fhPrimEtaAccPhi(0),\r
69fhPrimOmegaPt(0), fhPrimOmegaAccPt(0), fhPrimOmegaY(0), fhPrimOmegaAccY(0), \r
70fhPrimOmegaPhi(0), fhPrimOmegaAccPhi(0),\r
71fCalorimeter("")\r
72{\r
73 //Default Ctor\r
74 InitParameters();\r
75}\r
76\r
77//______________________________________________________________________________\r
78AliAnaNeutralMeson::AliAnaNeutralMeson(const AliAnaNeutralMeson & ex) : AliAnaPartCorrBaseClass(ex), \r
79fAnaPi0Eta(ex.fAnaPi0Eta), fAnaOmega(ex.fAnaOmega), fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv), fNAsy(ex.fNAsy),fNMod(ex.fNMod),\r
80fNHistos(ex.fNHistos), fPerPtBin(ex.fPerPtBin),\r
81fAsyCut(ex.fAsyCut), fModCut(ex.fModCut),\r
82fDist(ex.fDist), fAsy(ex.fAsy), fPt(ex.fPt), fMass(ex.fMass),\r
83fInvMassCut(ex.fInvMassCut), \r
84fNbinsAsy(ex.fNbinsAsy), fMinAsy(ex.fMinAsy), fMaxAsy(ex.fMaxAsy),\r
85fNbinsPt(ex.fNbinsPt), fMinPt(ex.fMinPt), fMaxPt(ex.fMaxPt),\r
86fNbinsM(ex.fNbinsM), fMinM(ex.fMinM), fMaxM(ex.fMaxM),\r
87fEventsListPi0Eta(ex.fEventsListPi0Eta), fEventsListOmega(ex.fEventsListOmega),\r
88fPi0Mass(ex.fPi0Mass), fPi0MassPeakWidthCut(ex.fPi0MassPeakWidthCut), fhNClusters(ex.fhNClusters),\r
89fhRecPhoton(ex.fhRecPhoton), fhRecPhotonEtaPhi(ex.fhRecPhotonEtaPhi),\r
90fReal2Gamma(ex.fReal2Gamma), fMix2Gamma(ex.fMix2Gamma),fRealOmega(ex.fRealOmega),\r
91fMixOmegaA(ex.fMixOmegaA),fMixOmegaB(ex.fMixOmegaB), fMixOmegaC(ex.fMixOmegaC),\r
92fRealTwoGammaAsyPtM(ex.fRealTwoGammaAsyPtM), fMixTwoGammaAsyPtM(ex.fMixTwoGammaAsyPtM),\r
93fRealPi0GammaAsyPtM(ex.fRealPi0GammaAsyPtM), fMixAPi0GammaAsyPtM(ex.fMixAPi0GammaAsyPtM),\r
94fMixBPi0GammaAsyPtM(ex.fMixBPi0GammaAsyPtM), fMixCPi0GammaAsyPtM(ex.fMixCPi0GammaAsyPtM),\r
95fhPrimPhotonPt(ex.fhPrimPhotonPt), fhPrimPhotonAccPt(ex.fhPrimPhotonAccPt), fhPrimPhotonY(ex.fhPrimPhotonY), \r
96fhPrimPhotonAccY(ex.fhPrimPhotonAccY), fhPrimPhotonPhi(ex.fhPrimPhotonPhi), \r
97fhPrimPhotonAccPhi(ex.fhPrimPhotonAccPhi),\r
98fhPrimPi0Pt(ex.fhPrimPi0Pt), fhPrimPi0AccPt(ex.fhPrimPi0AccPt), fhPrimPi0Y(ex.fhPrimPi0Y),\r
99fhPrimPi0AccY(ex.fhPrimPi0AccY), fhPrimPi0Phi(ex.fhPrimPi0Phi), fhPrimPi0AccPhi(ex.fhPrimPi0AccPhi), \r
100fhPrimEtaPt(ex.fhPrimEtaPt), fhPrimEtaAccPt(ex.fhPrimEtaAccPt), fhPrimEtaY(ex.fhPrimEtaY), \r
101fhPrimEtaAccY(ex.fhPrimEtaAccY), fhPrimEtaPhi(ex.fhPrimEtaPhi), fhPrimEtaAccPhi(ex.fhPrimEtaAccPhi),\r
102fhPrimOmegaPt(ex.fhPrimOmegaPt), fhPrimOmegaAccPt(ex.fhPrimOmegaAccPt), fhPrimOmegaY(ex.fhPrimOmegaY),\r
103fhPrimOmegaAccY(ex.fhPrimOmegaAccY), fhPrimOmegaPhi(ex.fhPrimOmegaPhi), fhPrimOmegaAccPhi(ex.fhPrimOmegaAccPhi),\r
104fCalorimeter("")\r
105{\r
106 // cpy ctor\r
107 //Do not need it\r
108}\r
109\r
110//______________________________________________________________________________\r
111AliAnaNeutralMeson & AliAnaNeutralMeson::operator = (const AliAnaNeutralMeson & ex)\r
112{\r
113 // assignment operator\r
114\r
115 if(this == &ex)return *this;\r
116 ((AliAnaPartCorrBaseClass *)this)->operator=(ex);\r
117 \r
118 fAnaPi0Eta = ex.fAnaPi0Eta; fAnaOmega = ex.fAnaOmega;\r
119 fNPID=ex.fNPID; fNmaxMixEv=ex.fNmaxMixEv; fNAsy=ex.fNAsy;fNMod=ex.fNMod;\r
120 fNHistos=ex.fNHistos; fPerPtBin=ex.fPerPtBin;\r
121 fAsyCut=ex.fAsyCut; fModCut=ex.fModCut;\r
122 fDist=ex.fDist; fAsy=ex.fAsy; fPt=ex.fPt; fMass=ex.fMass;\r
123 fInvMassCut=ex.fInvMassCut;\r
124 fNbinsAsy=ex.fNbinsAsy; fMinAsy=ex.fMinAsy; fMaxAsy=ex.fMaxAsy;\r
125 fNbinsPt=ex.fNbinsPt; fMinPt=ex.fMinPt; fMaxPt=ex.fMaxPt;\r
126 fNbinsM=ex.fNbinsM; fMinM=ex.fMinM; fMaxM=ex.fMaxM;\r
127 fEventsListPi0Eta=ex.fEventsListPi0Eta; fEventsListOmega=ex.fEventsListOmega;\r
128 fPi0Mass=ex.fPi0Mass; fPi0MassPeakWidthCut=ex.fPi0MassPeakWidthCut; fhNClusters =ex.fhNClusters;\r
129 fhRecPhoton=ex.fhRecPhoton; fhRecPhotonEtaPhi =ex.fhRecPhotonEtaPhi;\r
130 fReal2Gamma=ex.fReal2Gamma; fMix2Gamma=ex.fMix2Gamma;fRealOmega=ex.fRealOmega;\r
131 fMixOmegaA=ex.fMixOmegaA; fMixOmegaB=ex.fMixOmegaB; fMixOmegaC=ex.fMixOmegaC;\r
132 fRealTwoGammaAsyPtM=ex.fRealTwoGammaAsyPtM; fMixTwoGammaAsyPtM =ex.fMixTwoGammaAsyPtM;\r
133 fRealPi0GammaAsyPtM=ex.fRealPi0GammaAsyPtM; fMixAPi0GammaAsyPtM=ex.fMixAPi0GammaAsyPtM;\r
134 fMixBPi0GammaAsyPtM=ex.fMixBPi0GammaAsyPtM; fMixCPi0GammaAsyPtM=ex.fMixCPi0GammaAsyPtM;\r
135 fhPrimPhotonPt=ex.fhPrimPhotonPt; fhPrimPhotonAccPt=ex.fhPrimPhotonAccPt; fhPrimPhotonY=ex.fhPrimPhotonY;\r
136 fhPrimPhotonAccY=ex.fhPrimPhotonAccY; fhPrimPhotonPhi=ex.fhPrimPhotonPhi;\r
137 fhPrimPhotonAccPhi=ex.fhPrimPhotonAccPhi;\r
138 fhPrimPi0Pt=ex.fhPrimPi0Pt; fhPrimPi0AccPt=ex.fhPrimPi0AccPt; fhPrimPi0Y=ex.fhPrimPi0Y;\r
139 fhPrimPi0AccY=ex.fhPrimPi0AccY; fhPrimPi0Phi=ex.fhPrimPi0Phi; fhPrimPi0AccPhi=ex.fhPrimPi0AccPhi;\r
140 fhPrimEtaPt=ex.fhPrimEtaPt; fhPrimEtaAccPt=ex.fhPrimEtaAccPt; fhPrimEtaY=ex.fhPrimEtaY;\r
141 fhPrimEtaAccY=ex.fhPrimEtaAccY; fhPrimEtaPhi=ex.fhPrimEtaPhi; fhPrimEtaAccPhi=ex.fhPrimEtaAccPhi;\r
142 fhPrimOmegaPt=ex.fhPrimOmegaPt; fhPrimOmegaAccPt=ex.fhPrimOmegaAccPt; fhPrimOmegaY=ex.fhPrimOmegaY;\r
143 fhPrimOmegaAccY=ex.fhPrimOmegaAccY; fhPrimOmegaPhi=ex.fhPrimOmegaPhi; fhPrimOmegaAccPhi=ex.fhPrimOmegaAccPhi;\r
144 \r
145 return *this;\r
146 \r
147}\r
148\r
149//______________________________________________________________________________\r
150AliAnaNeutralMeson::~AliAnaNeutralMeson() {\r
151\r
152 //dtor\r
153 \r
154 if(fEventsListPi0Eta){\r
155 delete[] fEventsListPi0Eta;\r
156 fEventsListPi0Eta=0 ;\r
157 }\r
158 if(fEventsListOmega){\r
159 delete[] fEventsListOmega;\r
160 fEventsListOmega=0 ;\r
161 } \r
162\r
163#ifdef __PHOSGEO__\r
164 if(fPHOSGeo) delete fPHOSGeo ;\r
165#endif \r
166}\r
167\r
168//______________________________________________________________________________\r
169void AliAnaNeutralMeson::SetCalorimeter(TString det)\r
170{\r
171 //set the detector\r
172 //for PHOS, we consider three cases with different PHOS number \r
173 //for EMCAL, we use the default one, I don't know, needs to be confirmed\r
174 Int_t n=0;\r
175 Int_t * nmod=0;\r
176 if(det == "PHOS"){ \r
177 n = 3;\r
178 nmod= new Int_t [n];\r
179 nmod[0]=1; nmod[1]=3; nmod[2]=5;\r
180 }\r
181 \r
182 \r
183 if(det == "EMCAL") { //the default number of EMCAL modules\r
184 n = 1;\r
185 nmod= new Int_t [n]; \r
186 nmod[0]=0; \r
187 } \r
188 fCalorimeter = det;\r
189 fNMod = n;\r
190 fModCut = nmod;\r
191}\r
192\r
193//______________________________________________________________________________\r
194void AliAnaNeutralMeson::InitParameters()\r
195{\r
196//Init parameters when first called the analysis\r
197//Set default parameters\r
198 SetInputAODName("photons");\r
199 AddToHistogramsName("AnaNeutralMesons_");\r
200 \r
201 fInvMassCut = 1.;\r
202 \r
203 fAnaPi0Eta=kTRUE;\r
204 fAnaOmega=kTRUE;\r
205 \r
206 fNPID = 9; \r
207 fNAsy = 4;\r
208 fAsyCut = new Double_t[fNAsy] ;\r
209 fAsyCut[0]=0.7; fAsyCut[1]=0.8; fAsyCut[2]=0.9; fAsyCut[3]=1;\r
210 \r
211 fNHistos=20; // set the max pt range from 0 to 20 GeV\r
212 fPerPtBin=1.;//set the invariant mass spectrum in pt bin\r
213\r
214 fNmaxMixEv = 6; //event buffer size\r
215 fPi0Mass=0.135; //nominal pi0 mass in GeV/c^2\r
216 fPi0MassPeakWidthCut=0.015; //the Pi0MassCut should dependent on pt, here we assume it 0.015 \r
217 \r
218 fNbinsAsy =200; \r
219 fMinAsy =0;\r
220 fMaxAsy =1;\r
221 fNbinsPt =200;\r
222 fMinPt = 0;\r
223 fMaxPt=20;\r
224 fNbinsM = 200;\r
225 fMinM=0;\r
226 fMaxM=1.;\r
227}\r
228\r
229//______________________________________________________________________________\r
230void AliAnaNeutralMeson::Init()\r
231{ \r
232//Init some data members needed in analysis\r
233 \r
234 fEventsListPi0Eta = new TList ;\r
235 fEventsListOmega = new TList ;\r
236 \r
237#ifdef __PHOSGEO__\r
238 printf("PHOS geometry initialized!\n");\r
239 fPHOSGeo = new AliPHOSGeoUtils("PHOSgeo") ;\r
240#endif \r
241 \r
242}\r
243\r
244//______________________________________________________________________________\r
245TList * AliAnaNeutralMeson::GetCreateOutputObjects()\r
246{ \r
247 // Create histograms to be saved in output file and \r
248 // store them in fOutputContainer\r
249 \r
250 TList * outputContainer = new TList() ; \r
251 outputContainer->SetName(GetName()); \r
252 \r
253 char key[255] ;\r
254 char title[255] ;\r
255 const char * detector= fCalorimeter;\r
256\r
257 fhNClusters = new TH1I("fhNClusters","N clusters per event",100,0,100); //\r
258 outputContainer->Add(fhNClusters);\r
259\r
260 fhRecPhoton = new TH1D*[fNPID*fNMod];\r
261 fhRecPhotonEtaPhi= new TH2F*[fNPID*fNMod];\r
262 // plot the photon cluster pt, eta and phi distribution by using the PID and module cut \r
263 for(Int_t ipid=0;ipid<fNPID;ipid++){\r
264 for(Int_t jmod=0;jmod<fNMod;jmod++){\r
265 sprintf(key,"fhRecPhoton_%s_pid_%d_nmod_%d",detector,ipid,fModCut[jmod]);\r
266 sprintf(title,"Rec. photon with pid#%d in %d modules with %s",ipid, fModCut[jmod],detector );\r
267 fhRecPhoton[ipid*fNMod+jmod] = new TH1D(key,title,fNbinsPt,fMinPt, fMaxPt);\r
268 fhRecPhoton[ipid*fNMod+jmod]->GetXaxis()->SetTitle("pt (GeV/c)");\r
269 fhRecPhoton[ipid*fNMod+jmod]->GetYaxis()->SetTitle("Entries");\r
270 outputContainer->Add(fhRecPhoton[ipid*fNMod+jmod]);\r
271 \r
272 sprintf(key,"fhRecPhotonEtaPhi_%s_pid_%d_nmod_%d",detector, ipid,fModCut[jmod]);\r
273 sprintf(title,"Rec. photon eta vs phi with pid#%d in %d modules with %s",ipid, fModCut[jmod], detector);\r
274 fhRecPhotonEtaPhi[ipid*fNMod+jmod] = new TH2F(key,title,200,-1,1,200,0,7);\r
275 fhRecPhotonEtaPhi[ipid*fNMod+jmod]->GetXaxis()->SetTitle("#eta");\r
276 fhRecPhotonEtaPhi[ipid*fNMod+jmod]->GetYaxis()->SetTitle("#phi (rad.)");\r
277 outputContainer->Add(fhRecPhotonEtaPhi[ipid*fNMod+jmod]);\r
278 }\r
279 } \r
280\r
281\r
282 if(fAnaPi0Eta){\r
283 fReal2Gamma = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
284 fMix2Gamma = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
285 \r
286 fRealTwoGammaAsyPtM = new TH3F*[fNPID];\r
287 fMixTwoGammaAsyPtM = new TH3F*[fNPID];\r
288 for(Int_t ipid=0;ipid<fNPID;ipid++){\r
289 sprintf(key,"RealTwoGammaAsyPtM_%s_pid_%d",detector,ipid);\r
290 sprintf(title, "%s Real 2#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
291 fRealTwoGammaAsyPtM[ipid] = new TH3F(key, title,\r
292 fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
293 fRealTwoGammaAsyPtM[ipid]->GetYaxis()->SetTitle("2#gamma pt");\r
294 fRealTwoGammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#gamma1}-E_{#gamma2}|/(E_{#gamma1}+E_{#gamma2})");\r
295 fRealTwoGammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{2#gamma}");\r
296 outputContainer->Add(fRealTwoGammaAsyPtM[ipid]);\r
297 \r
298 sprintf(key,"MixTwoGammaAsyPtM_%s_pid_%d",detector,ipid);\r
299 sprintf(title, "%s Mix 2#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
300 fMixTwoGammaAsyPtM[ipid] = new TH3F(key, title,\r
301 fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
302 fMixTwoGammaAsyPtM[ipid]->GetYaxis()->SetTitle("2#gamma pt");\r
303 fMixTwoGammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#gamma1}-E_{#gamma2}|/(E_{#gamma1}+E_{#gamma2})");\r
304 fMixTwoGammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{2#gamma}");\r
305 outputContainer->Add(fMixTwoGammaAsyPtM[ipid]); \r
306 \r
307 for(Int_t jmod=0; jmod<fNMod; jmod++){\r
308 for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
309 for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
310 Double_t pt1=mhist*fPerPtBin;\r
311 Double_t pt2=(mhist+1)*fPerPtBin;\r
312 Int_t index = ipid*fNMod*fNAsy*fNHistos+jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
313 \r
314 sprintf(key,"fReal2Gamma_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector, ipid,fModCut[jmod],kasy,mhist);\r
315 sprintf(title,"real 2#gamma IVM with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
316 ipid,fAsyCut[kasy], fModCut[jmod],detector, pt1, pt2);\r
317 fReal2Gamma[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
318 fReal2Gamma[index]->GetYaxis()->SetTitle("dN/dM_{#gamma#gamma}");\r
319 fReal2Gamma[index]->GetXaxis()->SetTitle("M_{#gamma#gamma} (GeV/c^{2})");\r
320 outputContainer->Add(fReal2Gamma[index]);\r
321 \r
322 sprintf(key,"fMix2Gamma_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector,ipid,fModCut[jmod],kasy,mhist);\r
323 sprintf(title,"Mix 2#gamma IVM with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
324 ipid,fAsyCut[kasy],fModCut[jmod],detector,pt1, pt2);\r
325 fMix2Gamma[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
326 fMix2Gamma[index]->GetYaxis()->SetTitle("dN/dM_{#gamma#gamma}");\r
327 fMix2Gamma[index]->GetXaxis()->SetTitle("M_{#gamma#gamma} (GeV/c^{2})");\r
328 outputContainer->Add(fMix2Gamma[index]);\r
329 } \r
330 }\r
331 }\r
332 }\r
333 }\r
334 \r
335 if(fAnaOmega){\r
336 fRealOmega = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
337 fMixOmegaA = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
338 fMixOmegaB = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
339 fMixOmegaC = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
340 fRealPi0GammaAsyPtM=new TH3F*[fNPID];\r
341 fMixAPi0GammaAsyPtM=new TH3F*[fNPID];\r
342 fMixBPi0GammaAsyPtM=new TH3F*[fNPID];\r
343 fMixCPi0GammaAsyPtM=new TH3F*[fNPID];\r
344 \r
345 for(Int_t ipid=0;ipid<fNPID;ipid++){\r
346 sprintf(key,"RealPi0GammaAsyPtM_%s_pid_%d",detector,ipid);\r
347 sprintf(title,"%s Real #omega->#pi^{0}#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
348 fRealPi0GammaAsyPtM[ipid] = new TH3F(key, title,\r
349 fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
350 fRealPi0GammaAsyPtM[ipid]->GetYaxis()->SetTitle("#pi^{0}#gamma pt");\r
351 fRealPi0GammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#pi^{0}}-E_{#gamma}|/(E_{#pi^{0}}+E_{#gamma})");\r
352 fRealPi0GammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{#pi^{0}#gamma}");\r
353 outputContainer->Add(fRealPi0GammaAsyPtM[ipid]);\r
354 \r
355 sprintf(key,"MixAPi0GammaAsyPtM_%s_pid_%d",detector, ipid);\r
356 sprintf(title,"%s MixA #omega->#pi^{0}#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
357 fMixAPi0GammaAsyPtM[ipid] = new TH3F(key,title,\r
358 fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
359 fMixAPi0GammaAsyPtM[ipid]->GetYaxis()->SetTitle("#pi^{0}#gamma pt");\r
360 fMixAPi0GammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#pi^{0}}-E_{#gamma}|/(E_{#pi^{0}}+E_{#gamma})");\r
361 fMixAPi0GammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{#pi^{0}#gamma}");\r
362 outputContainer->Add(fMixAPi0GammaAsyPtM[ipid]);\r
363\r
364 sprintf(key,"MixBPi0GammaAsyPtM_%s_pid_%d",detector,ipid);\r
365 sprintf(title,"%s MixB #omega->#pi^{0}#gamma asy:pt:ivm with pid#%d",detector,ipid); \r
366 fMixBPi0GammaAsyPtM[ipid] = new TH3F(key, title,\r
367 fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
368 fMixBPi0GammaAsyPtM[ipid]->GetYaxis()->SetTitle("#pi^{0}#gamma pt");\r
369 fMixBPi0GammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#pi^{0}}-E_{#gamma}|/(E_{#pi^{0}}+E_{#gamma})");\r
370 fMixBPi0GammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{#pi^{0}#gamma}");\r
371 outputContainer->Add(fMixBPi0GammaAsyPtM[ipid]);\r
372 \r
373 sprintf(key,"MixCPi0GammaAsyPtM_%s_pid_%d",detector,ipid);\r
374 sprintf(title,"%s MixC #omega->#pi^{0}#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
375 fMixCPi0GammaAsyPtM[ipid] = new TH3F(key, title,\r
376 fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
377 fMixCPi0GammaAsyPtM[ipid]->GetYaxis()->SetTitle("#pi^{0}#gamma pt");\r
378 fMixCPi0GammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#pi^{0}}-E_{#gamma}|/(E_{#pi^{0}}+E_{#gamma})");\r
379 fMixCPi0GammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{#pi^{0}#gamma}");\r
380 outputContainer->Add(fMixCPi0GammaAsyPtM[ipid]);\r
381\r
382 for(Int_t jmod=0; jmod<fNMod; jmod++){\r
383 for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
384 for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
385 Double_t pt1=mhist*fPerPtBin;\r
386 Double_t pt2=(mhist+1)*fPerPtBin;\r
387 Int_t index = ipid*fNMod*fNAsy*fNHistos+jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
388 sprintf(key,"fRealOmega_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector,ipid,fModCut[jmod],kasy,mhist);\r
389 sprintf(title,"Real #omega with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
390 ipid,fAsyCut[kasy],fModCut[jmod], detector, pt1, pt2);\r
391 fRealOmega[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
392 fRealOmega[index]->GetYaxis()->SetTitle("dN/dM_{#pi^{0}#gamma} "); \r
393 fRealOmega[index]->GetXaxis()->SetTitle("M_{#pi^{0}#gamma} (GeV/c^{2})");\r
394 outputContainer->Add(fRealOmega[index]);\r
395 \r
396 sprintf(key,"fMixOmegaA_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector, ipid,fModCut[jmod],kasy,mhist);\r
397 sprintf(title,"Mix #omega A with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
398 ipid,fAsyCut[kasy],fModCut[jmod],detector,pt1, pt2);\r
399 fMixOmegaA[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
400 fMixOmegaA[index]->GetYaxis()->SetTitle("dN/dM_{#pi^{0}#gamma}");\r
401 fMixOmegaA[index]->GetXaxis()->SetTitle("M_{#pi^{0}#gamma} (GeV/c^{2})");\r
402 outputContainer->Add(fMixOmegaA[index]);\r
403 \r
404 sprintf(key,"fMixOmegaB_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector,ipid,fModCut[jmod],kasy,mhist);\r
405 sprintf(title,"Mix #omega B with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
406 ipid,fAsyCut[kasy],fModCut[jmod],detector,pt1, pt2);\r
407 fMixOmegaB[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
408 fMixOmegaB[index]->GetYaxis()->SetTitle("dN/dM_{#pi^{0}#gamma}");\r
409 fMixOmegaB[index]->GetXaxis()->SetTitle("M_{#pi^{0}#gamma} (GeV/c^{2})");\r
410 outputContainer->Add(fMixOmegaB[index]);\r
411 \r
412 sprintf(key,"fMixOmegaC_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector,ipid,fModCut[jmod],kasy,mhist);\r
413 sprintf(title,"Mix #omega C with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
414 ipid,fAsyCut[kasy],fModCut[jmod],detector, pt1, pt2);\r
415 fMixOmegaC[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
416 fMixOmegaC[index]->GetYaxis()->SetTitle("dN/dM_{#pi^{0}#gamma}");\r
417 fMixOmegaC[index]->GetXaxis()->SetTitle("M_{#pi^{0}#gamma} (GeV/c^{2})");\r
418 outputContainer->Add(fMixOmegaC[index]);\r
419 }\r
420 }\r
421 }\r
422 }\r
423 }\r
424\r
425 //Histograms filled only if MC data is requested \r
426 //pt distribution\r
427 //currently only PHOS included \r
428 if(fCalorimeter=="PHOS" && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){\r
429 fhPrimPhotonPt = new TH1D("hPrimPhotonPt","Primary phton pt at |y|=1",fNbinsPt,fMinPt, fMaxPt) ;\r
430 outputContainer->Add(fhPrimPhotonPt);\r
431\r
432 fhPrimPhotonAccPt = new TH1D("hPrimPhotonAccPt","Primary photon pt in acceptance",fNbinsPt,fMinPt, fMaxPt);\r
433 outputContainer->Add(fhPrimPhotonAccPt);\r
434 \r
435 fhPrimPi0Pt = new TH1D("hPrimPi0Pt","Primary pi0 pt at |y|=1",fNbinsPt,fMinPt, fMaxPt) ;\r
436 outputContainer->Add(fhPrimPi0Pt);\r
437\r
438 fhPrimPi0AccPt = new TH1D("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",fNbinsPt,fMinPt, fMaxPt);\r
439 outputContainer->Add(fhPrimPi0AccPt);\r
440 \r
441 fhPrimEtaPt = new TH1D("hPrimEtaPt","Primary eta pt at |y|=1",fNbinsPt,fMinPt, fMaxPt) ;\r
442 outputContainer->Add(fhPrimEtaPt);\r
443\r
444 fhPrimEtaAccPt = new TH1D("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",fNbinsPt,fMinPt, fMaxPt) ;\r
445 outputContainer->Add(fhPrimEtaAccPt);\r
446 \r
447 fhPrimOmegaPt = new TH1D("hPrimOmegaPt","Primary Omega pt at |y|=1",fNbinsPt,fMinPt, fMaxPt) ;\r
448 outputContainer->Add(fhPrimOmegaPt);\r
449\r
450 fhPrimOmegaAccPt = new TH1D("hPrimOmegaAccPt","Primary Omega pt with omega->pi0+gamma in acceptance",fNbinsPt,fMinPt, fMaxPt);\r
451 outputContainer->Add(fhPrimOmegaAccPt);\r
452 \r
453 //for photon\r
454 fhPrimPhotonY = new TH1D("hPrimPhotonY","Rapidity of primary pi0",100,-5.,5.) ;\r
455 outputContainer->Add(fhPrimPhotonY) ;\r
456 \r
457 fhPrimPhotonAccY = new TH1D("hPrimPhotonAccY","Rapidity of primary pi0",100,-5.,5.) ;\r
458 outputContainer->Add(fhPrimPhotonAccY) ;\r
459 \r
460 fhPrimPhotonPhi = new TH1D("hPrimPhotonPhi","Azimithal of primary photon",180,0.,360.) ;\r
461 outputContainer->Add(fhPrimPhotonPhi) ;\r
462 \r
463 fhPrimPhotonAccPhi = new TH1D("hPrimPhotonAccPhi","Azimithal of primary photon in Acceptance",180,-0.,360.) ;\r
464 outputContainer->Add(fhPrimPhotonAccPhi) ;\r
465 \r
466 //for pi0\r
467 fhPrimPi0Y = new TH1D("hPrimPi0Y","Rapidity of primary pi0",100,-5.,5.) ;\r
468 outputContainer->Add(fhPrimPi0Y) ;\r
469 \r
470 fhPrimPi0AccY = new TH1D("hPrimPi0AccY","Rapidity of primary pi0",100,-5.,5.) ;\r
471 outputContainer->Add(fhPrimPi0AccY) ;\r
472 \r
473 fhPrimPi0Phi = new TH1D("hPrimPi0Phi","Azimithal of primary pi0",180,0.,360.) ;\r
474 outputContainer->Add(fhPrimPi0Phi) ;\r
475 \r
476 fhPrimPi0AccPhi = new TH1D("hPrimPi0AccPhi","Azimithal of primary pi0 with accepted daughters",180,-0.,360.) ;\r
477 outputContainer->Add(fhPrimPi0AccPhi) ;\r
478 \r
479 //for eta\r
480 fhPrimEtaY = new TH1D("hPrimEtaY","Rapidity of primary Eta",100,-5.,5.) ;\r
481 outputContainer->Add(fhPrimEtaY) ;\r
482 \r
483 fhPrimEtaAccY = new TH1D("hPrimEtaAccY","Rapidity of primary eta",100,-5.,5.) ;\r
484 outputContainer->Add(fhPrimEtaAccY) ;\r
485 \r
486 fhPrimEtaPhi = new TH1D("hPrimEtaPhi","Azimithal of primary eta",180,0.,360.) ;\r
487 outputContainer->Add(fhPrimEtaPhi) ;\r
488 \r
489 fhPrimEtaAccPhi = new TH1D("hPrimEtaAccPhi","Azimithal of primary eta with accepted daughters",180,-0.,360.) ;\r
490 outputContainer->Add(fhPrimEtaAccPhi) ;\r
491 \r
492 //for omega\r
493 fhPrimOmegaY = new TH1D("hPrimOmegaY","Rapidity of primary Omega",100,-5.,5.) ;\r
494 outputContainer->Add(fhPrimOmegaY) ;\r
495 \r
496 fhPrimOmegaAccY = new TH1D("hPrimOmegaAccY","Rapidity of primary omega->pi0+gamma",100,-5.,5.) ;\r
497 outputContainer->Add(fhPrimOmegaAccY) ;\r
498 \r
499 fhPrimOmegaPhi = new TH1D("hPrimOmegaPhi","Azimithal of primary omega",180,0.,360.) ;\r
500 outputContainer->Add(fhPrimOmegaPhi) ;\r
501 \r
502 fhPrimOmegaAccPhi = new TH1D("hPrimOmegaAccPhi","Azimithal of primary omega->pi0+gamma with accepted daughters",180,-0.,360.) ;\r
503 outputContainer->Add(fhPrimOmegaAccPhi);\r
504 }\r
505 return outputContainer;\r
506}\r
507\r
508//______________________________________________________________________________\r
509void AliAnaNeutralMeson::Print(const Option_t * /*opt*/) const\r
510{\r
511 //Print some relevant parameters set for the analysis\r
512 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
513 AliAnaPartCorrBaseClass::Print(" ");\r
514 printf("performing the analysis of 2gamma IVM %d\n", fAnaPi0Eta); //IVM: Invariant Mass\r
515 printf("performing the analysis of pi0+gamma->3gamma IVM %d\n", fAnaOmega);\r
516 printf("Number of Events to be mixed: %d \n",fNmaxMixEv) ;\r
517 printf("Number of different PID used: %d \n",fNPID) ;\r
518 printf("Number of different Asy cut: %d\n",fNAsy);\r
519 printf("Asy bins: %d Min: %2.1f Max: %2.1f\n", fNbinsAsy, fMinAsy, fMaxAsy);\r
520 printf("Pt bins: %d Min: %2.1f Max: %2.1f GeV/c\n", fNbinsPt, fMinPt, fMaxPt);\r
521 printf("IVM bins: %d Min: %2.1f Max: %2.1f GeV/c^2\n", fNbinsM, fMinM, fMaxM);\r
522 printf("produce the %d histograms per %2.1fGeV/c pt bin \n", fNHistos,fPerPtBin);\r
523 printf("Vertec, centrality Cuts, RP are not used at present \n") ;\r
524 printf("------------------------------------------------------\n") ;\r
525\r
526} \r
527\r
528//______________________________________________________________________________\r
529void AliAnaNeutralMeson::MakeAnalysisFillHistograms() \r
530{\r
531 //process event from AOD brach\r
532 //extract pi0, eta and omega analysis\r
533 Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;\r
534 if(IsBadRun(iRun)) return ; \r
535\r
536 //vertex cut not used yet\r
537 //centrality not used yet\r
538 //reaction plane not used yet\r
539\r
540 TClonesArray * array1 = (TClonesArray*)GetInputAODBranch();\r
541 Int_t nPhot = array1 ->GetEntries();\r
542 fhNClusters->Fill(nPhot);\r
543\r
544 //fill photon clusters\r
545 for(Int_t i=0;i<nPhot;i++){\r
546 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i)) ;\r
547 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
548 for(Int_t ipid=0;ipid<fNPID;ipid++){\r
549 for(Int_t jmod=0;jmod<fNMod;jmod++){\r
550 if(IsPhotonSelected(p1, ipid, fModCut[jmod])){\r
551 fhRecPhoton[ipid*fNMod+jmod]->Fill(photon1.Pt());\r
552 Double_t phi = photon1.Phi() ;\r
553 if(phi<0) phi = photon1.Phi()+2*TMath::Pi();\r
554 fhRecPhotonEtaPhi[ipid*fNMod+jmod]->Fill(photon1.Eta(),phi);\r
555 } \r
556 }\r
557 }\r
558 }\r
559\r
560/////////////////////////////////////////////\r
561//real events\r
562 if(nPhot>=2 && fAnaPi0Eta){\r
563// printf("real nPhot: %d \n",nPhot);\r
564 \r
565 RealPi0Eta(array1); //real events for pi0 and eta\r
566 \r
567 TList * evMixListPi0Eta=fEventsListPi0Eta ; //with cluster larger than 1\r
568 Int_t nMixed = evMixListPi0Eta->GetSize() ;\r
569 ////////////////////////////////////////////\r
570 //Mixing events for pi0 and eta to reconstruct the background\r
571 for(Int_t i1=0; i1<nMixed; i1++){\r
572 TClonesArray* array2= (TClonesArray*) (evMixListPi0Eta->At(i1));\r
573 MixPi0Eta(array1, array2);\r
574 }\r
575 ////////////////////////////////////////////\r
576 \r
577 //event buffer for pi0 and eta\r
578 TClonesArray *currentEvent = new TClonesArray(*array1);\r
579 \r
580 if(currentEvent->GetEntriesFast()>=2){\r
581 evMixListPi0Eta->AddFirst(currentEvent) ;\r
582 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer\r
583 if(evMixListPi0Eta->GetSize()>=fNmaxMixEv) {\r
584 TClonesArray * tmp = (TClonesArray*) (evMixListPi0Eta->Last()) ;\r
585 evMixListPi0Eta->RemoveLast() ;\r
586 delete tmp ;\r
587 }\r
588 }\r
589 else{ //empty event\r
590 delete currentEvent ;\r
591 currentEvent=0 ;\r
592 }\r
593 \r
594 }\r
595 \r
596 if(nPhot>=3 && fAnaOmega) {\r
597 RealOmega(array1); //real events for omega\r
598 TList * evMixListOmega = fEventsListOmega ; //with cluster larger than 2\r
599 Int_t nMixedOmega = evMixListOmega->GetSize() ;\r
600 \r
601 ////////////////////////////////////////////\r
602 //mix for omega\r
603 for(Int_t i1=0; i1<nMixedOmega; i1++){\r
604 TClonesArray* array2= (TClonesArray*) (evMixListOmega->At(i1));\r
605 MixOmegaAB(array1, array2);\r
606 //third event\r
607 for(Int_t i2=i1+1;i2<nMixedOmega;i2++){\r
608 TClonesArray * array3 = (TClonesArray*)(evMixListOmega->At(i2));\r
609// printf("omega nPhot: %d nPhot2: %d nPhot3: %d \n", nPhot,nPhot2,nPhot3);\r
610 MixOmegaC(array1,array2,array3);\r
611 }\r
612 }\r
613 \r
614 //fill event buffer with 2 more clusters events for omega extraction\r
615 TClonesArray *currentEventOmega = new TClonesArray(*array1);\r
616 if(currentEventOmega->GetEntriesFast()>=3){\r
617 evMixListOmega->AddFirst(currentEventOmega) ;\r
618 currentEventOmega=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer\r
619 if(evMixListOmega->GetSize()>=fNmaxMixEv) {\r
620 TClonesArray * tmp = (TClonesArray*) (evMixListOmega->Last()) ;\r
621 evMixListOmega->RemoveLast() ;\r
622 delete tmp ;\r
623 }\r
624 } \r
625 else{ //empty event\r
626 delete currentEventOmega ;\r
627 currentEventOmega=0 ;\r
628 } \r
629 \r
630 } \r
631\r
632 if(fCalorimeter=="PHOS" && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){\r
633 AliStack * stack = GetMCStack();\r
634 //photon acceptance\r
635 PhotonAcceptance(stack);\r
636 //pi0 and eta acceptance\r
637 Pi0EtaAcceptance(stack);\r
638 //omega acceptance \r
639 OmegaAcceptance(stack);\r
640 }\r
641\r
642}\r
643\r
644//______________________________________________________________________________\r
645Bool_t AliAnaNeutralMeson::IsPhotonSelected(AliAODPWG4Particle *p1, Int_t ipid, Int_t nmod)\r
646{\r
647 //select the photon to be analyzed based on pid and which module \r
648\r
649 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
650\r
651 if(fCalorimeter == "PHOS"){\r
652 Int_t mod = 20 ;\r
653#ifdef __PHOSGEO__\r
654 Double_t vtx[3]={0,0,0}; \r
655// if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vtx);\r
656 Double_t x = 0 ,z = 0 ;\r
657 fPHOSGeo->ImpactOnEmc(vtx,p1->Theta(),p1->Phi(), mod,z,x) ;\r
658#endif\r
659 if(mod<=nmod && p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) {return kTRUE;}\r
660 }\r
661 \r
662 else if(fCalorimeter == "EMCAL" && p1->IsPIDOK(ipid,AliCaloPID::kPhoton) ) {return kTRUE;} //default EMCAL module\r
663 return kFALSE;\r
664}\r
665\r
666//______________________________________________________________________________\r
667void AliAnaNeutralMeson::Get2GammaAsyDistPtM(AliAODPWG4Particle *p1, AliAODPWG4Particle *p2,\r
668 Double_t &asy, Double_t &dist, Double_t &pt, Double_t &mass)\r
669{\r
670 //get the asy, dist, pair pt and pair inv mass from the 2gammas\r
671 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
672 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
673 if(photon1.Pt()>0 && photon2.Pt()>0){\r
674 asy=TMath::Abs(photon1.E()-photon2.E())/(photon1.E()+photon2.E());\r
675 Double_t phi1 = photon1.Phi();\r
676 Double_t phi2 = photon2.Phi();\r
677 Double_t eta1 = photon1.Eta();\r
678 Double_t eta2 = photon2.Eta();\r
679 dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));\r
680 pt=(photon1+photon2).Pt();\r
681 mass=(photon1+photon2).M();\r
682// printf("asy: %2.3f dist: %2.3f pt: %2.3f m: %2.3f\n",asy,dist,pt,mass);\r
683 }\r
684 else { asy=0; dist =0; pt =0; mass=0;}\r
685}\r
686\r
687//______________________________________________________________________________\r
688void AliAnaNeutralMeson::GetAsyDistPtM(TLorentzVector photon1, TLorentzVector photon2,\r
689 Double_t &asy, Double_t &dist, Double_t &pt, Double_t &mass)\r
690{\r
691 //get the asy, dist, pair pt and pair inv mass from the 2gammas \r
692 if(photon1.Pt()>0 && photon2.Pt()>0){\r
693 asy=TMath::Abs(photon1.E()-photon2.E())/(photon1.E()+photon2.E());\r
694 Double_t phi1 = photon1.Phi();\r
695 Double_t phi2 = photon2.Phi();\r
696 Double_t eta1 = photon1.Eta();\r
697 Double_t eta2 = photon2.Eta();\r
698 dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));\r
699 pt=(photon1+photon2).Pt();\r
700 mass=(photon1+photon2).M();\r
701// printf("asy: %2.3f dist: %2.3f pt: %2.3f m: %2.3f\n",asy,dist,pt,mass);\r
702 }\r
703 else { asy=0; dist =0; pt =0; mass=0;}\r
704}\r
705\r
706//______________________________________________________________________________ \r
707void AliAnaNeutralMeson::RealPi0Eta(TClonesArray * array)\r
708{\r
709 //get the 2gamma invariant mass distribution from one event\r
710 Int_t nclusters = array->GetEntries();\r
711 for(Int_t ie=0;ie<nclusters;ie++){\r
712 for(Int_t je=ie+1;je<nclusters;je++){\r
713 AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array->At(ie);\r
714 AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array->At(je);\r
715 //Asy vs. pt\r
716 Get2GammaAsyDistPtM(p1,p2,fAsy,fDist,fPt,fMass);\r
717 for(Int_t ipid=0;ipid<fNPID;ipid++){\r
718\r
719 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ \r
720 fRealTwoGammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module with pid\r
721 }\r
722\r
723 for(Int_t jmod=0;jmod<fNMod;jmod++){\r
724 for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
725 if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && IsPhotonSelected(p2,ipid,fModCut[jmod])){ \r
726 Get2GammaAsyDistPtM(p1,p2,fAsy,fDist,fPt,fMass);\r
727 //filling the histo\r
728 for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
729 Double_t pt1 =mhist *fPerPtBin;\r
730 Double_t pt2 =(mhist+1)*fPerPtBin;\r
731 Int_t index = ipid*fNMod*fNAsy*fNHistos+jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
732 if(fMass<fInvMassCut && fAsy<fAsyCut[kasy] \r
733 && fPt<pt2 && fPt>pt1 &&fAsy && fDist && fPt && fMass)\r
734 fReal2Gamma[index]->Fill(fMass);\r
735 }//hist\r
736 } //photon selection \r
737 }//asy\r
738 }//mod\r
739 }//pid\r
740 }//je\r
741 }//ie\r
742 \r
743}\r
744\r
745//______________________________________________________________________________ \r
746void AliAnaNeutralMeson::MixPi0Eta(TClonesArray * array1, TClonesArray * array2)\r
747{\r
748 //mixing events for pi0 and eta invariant mass analysis\r
749 // printf("mixed events for pi0 and eta invariant mass analysis!\n");\r
750 Int_t nclusters1 = array1->GetEntries();\r
751 Int_t nclusters2 = array2->GetEntries();\r
752 for(Int_t ie=0;ie<nclusters1;ie++){\r
753 AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array1->At(ie);\r
754 for(Int_t je=0;je<nclusters2;je++){\r
755 AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array2->At(je);\r
756 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
757 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
758 Get2GammaAsyDistPtM(p1,p2, fAsy,fDist,fPt,fMass);\r
759\r
760 for(Int_t ipid=0;ipid<fNPID;ipid++){\r
761 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){\r
762 fMixTwoGammaAsyPtM[ipid]->Fill(fAsy,fPt,fMass); //default module with pid\r
763 }\r
764 for(Int_t jmod=0;jmod<fNMod;jmod++){ \r
765 for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
766 if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && IsPhotonSelected(p2,ipid,fModCut[jmod])){\r
767 for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
768 Double_t pt1 =mhist *fPerPtBin;\r
769 Double_t pt2 =(mhist+1)*fPerPtBin;\r
770 Int_t index = ipid*fNMod*fNAsy*fNHistos+jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
771 if(fMass<fInvMassCut && fAsy<fAsyCut[kasy] \r
772 && fPt<pt2 && fPt>pt1 &&fAsy&& fDist && fPt && fMass)\r
773 fMix2Gamma[index]->Fill(fMass);\r
774 } //hist\r
775 } //photon selected\r
776 } //asy\r
777 } //mod\r
778 }//pid\r
779 } //cluster2\r
780 } //cluster1 \r
781} \r
782\r
783//______________________________________________________________________________ \r
784void AliAnaNeutralMeson::RealOmega(TClonesArray * array)\r
785{\r
786 //omega invariant mass extraction\r
787 //\r
788 // 1. select the pi0 candidate from 2gamma \r
789 // 2. combine the third photon with the pi0 candidate to get the omega invariant mass\r
790 //\r
791 Int_t nclusters = array->GetEntries();\r
792 for(Int_t ie=0;ie<nclusters;ie++){//ptc1\r
793 for(Int_t je=ie+1;je<nclusters;je++){ //ptc2\r
794 AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array->At(ie);\r
795 AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array->At(je);\r
796 //pid, mod and asy cut\r
797 Double_t f2gammaAsy;\r
798 Get2GammaAsyDistPtM(p1,p2,f2gammaAsy,fDist,fPt,fMass);\r
799\r
800 if(TMath::Abs(fMass-fPi0Mass)<fPi0MassPeakWidthCut){//select the pi0 candidate\r
801 for(Int_t ke=0;ke<nclusters;ke++){ //ptc3\r
802 AliAODPWG4Particle * p3 = (AliAODPWG4Particle *)array->At(ke); \r
803 if(p3 != p2 && p3 != p1 ){ ////p3!=p1 p3!=p2\r
804 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
805 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
806 Double_t px12 =(photon1+photon2).Px();\r
807 Double_t py12 =(photon1+photon2).Py();\r
808 Double_t pz12 =(photon1+photon2).Pz();\r
809 Double_t e12 =(photon1+photon2).E();\r
810 TLorentzVector photon12(px12, py12,pz12,e12);\r
811 \r
812 TLorentzVector photon3(p3->Px(),p3->Py(),p3->Pz(),p3->E());\r
813 GetAsyDistPtM(photon12, photon3,fAsy,fDist,fPt,fMass);\r
814 for(Int_t ipid=0;ipid<fNPID;ipid++){ //pid\r
815 if(p1->IsPIDOK(ipid,AliCaloPID::kPhoton) && \r
816 p2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
817 p3->IsPIDOK(ipid,AliCaloPID::kPhoton)){\r
818 fRealPi0GammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module, without pid\r
819 }\r
820\r
821 for(Int_t jmod=0;jmod<fNMod;jmod++){ //mod\r
822 for(Int_t kasy=0;kasy<fNAsy;kasy++){ //asy\r
823 if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && \r
824 IsPhotonSelected(p2,ipid,fModCut[jmod]) &&\r
825 IsPhotonSelected(p3,ipid,fModCut[jmod])){ //photon selection\r
826 \r
827 //filling the histo\r
828 for(Int_t mhist=0;mhist<fNHistos;mhist++){ //hist\r
829 Double_t pt1 =mhist *fPerPtBin;\r
830 Double_t pt2 =(mhist+1)*fPerPtBin;\r
831 Int_t index = ipid*fNMod*fNAsy*fNHistos+\r
832 jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
833 if(fMass<fInvMassCut && f2gammaAsy<fAsyCut[kasy] && \r
834 fPt<pt2 && fPt>pt1&& fAsy<fAsyCut[kasy] &&fAsy&&fDist&&fPt&&fMass)\r
835 fRealOmega[index]->Fill(fMass);\r
836 }//hist\r
837 } //photon selection\r
838 } //asy\r
839 }//mod\r
840 }//pid\r
841 } //p3!=p1 p3!=p2\r
842 } //ptc3\r
843 } //pi0 candidate\r
844 }//je ptc2 \r
845 } //ie ptc1\r
846}\r
847\r
848\r
849//______________________________________________________________________________\r
850void AliAnaNeutralMeson::MixOmegaAB(TClonesArray * array1, TClonesArray * array2)\r
851{\r
852 //omega background\r
853 //three omega background, we classify it A, B and C\r
854 // --A\r
855 // (r1_event1+r2_event1)+r3_event2\r
856 //\r
857 // --B\r
858 // (r1_event1+r2_event2)+r3_event2\r
859 // \r
860 // --C\r
861 // (r1_event1+r2_event2)+r3_event3\r
862 //\r
863 Int_t nclusters1 = array1->GetEntries();\r
864 Int_t nclusters2 = array2->GetEntries();\r
865 for(Int_t ie=0;ie<nclusters1;ie++){\r
866 //mix Omega A\r
867 for(Int_t je=ie+1;je<nclusters1;je++){ //.............\r
868 AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array1->At(ie);\r
869 AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array1->At(je);\r
870 Double_t f2gammaAsy ;\r
871 Get2GammaAsyDistPtM(p1,p2,f2gammaAsy,fDist,fPt,fMass);\r
872 if(TMath::Abs(fMass-fPi0Mass)<fPi0MassPeakWidthCut ){ //pi0 candidate \r
873 for(Int_t ke=0;ke<nclusters2;ke++){ //third photon\r
874 AliAODPWG4Particle * p3 = (AliAODPWG4Particle *)array2->At(ke);\r
875 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
876 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
877 Double_t px12 =(photon1+photon2).Px();\r
878 Double_t py12 =(photon1+photon2).Py();\r
879 Double_t pz12 =(photon1+photon2).Pz();\r
880 Double_t e12 =(photon1+photon2).E();\r
881 TLorentzVector photon12(px12, py12,pz12,e12);\r
882 TLorentzVector photon3(p3->Px(),p3->Py(),p3->Pz(),p3->E());\r
883 GetAsyDistPtM(photon12, photon3,fAsy,fDist,fPt,fMass);\r
884 //pid, mod and asy cut\r
885 for(Int_t ipid=0;ipid<fNPID;ipid++){ //pid\r
886 if(p1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
887 p2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
888 p3->IsPIDOK(ipid,AliCaloPID::kPhoton)){\r
889 fMixAPi0GammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module, with pid\r
890 }\r
891 for(Int_t jmod=0;jmod<fNMod;jmod++){ //mod\r
892 for(Int_t kasy=0;kasy<fNAsy;kasy++){ //asy\r
893 if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && \r
894 IsPhotonSelected(p2,ipid,fModCut[jmod]) &&\r
895 IsPhotonSelected(p3,ipid,fModCut[jmod])){ //photon selection\r
896 //fill the hist\r
897 for(Int_t mhist=0;mhist<fNHistos;mhist++){ //hist\r
898 Double_t pt1 =mhist *fPerPtBin;\r
899 Double_t pt2 =(mhist+1)*fPerPtBin;\r
900 Int_t index = ipid*fNMod*fNAsy*fNHistos+ jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
901 if(fMass<fInvMassCut && f2gammaAsy<fAsyCut[kasy] &&\r
902 fPt<pt2 && fPt>pt1&& fAsy<fAsyCut[kasy] &&fAsy&&fDist&&fPt&&fMass)\r
903 fMixOmegaA[index]->Fill(fMass);\r
904 }//hist\r
905 }//photon selection\r
906 } //asy\r
907 }//mod\r
908 } //pid\r
909 } //ptc3\r
910 } //pi0 candidate\r
911 } //mixA je ptc2\r
912\r
913 //mix omega B\r
914 for(Int_t je=0;je<nclusters2;je++){ //mixb...............\r
915 AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array1->At(ie);\r
916 AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array2->At(je);\r
917 Double_t f2gammaAsy;\r
918 Get2GammaAsyDistPtM(p1,p2,f2gammaAsy,fDist,fPt,fMass);\r
919 if(TMath::Abs(fMass-fPi0Mass)<fPi0MassPeakWidthCut){\r
920 for(Int_t ke=0;ke<nclusters2;ke++){ //third photon\r
921 AliAODPWG4Particle * p3 = (AliAODPWG4Particle *)array2->At(ke);\r
922 if(p3 != p2){ //p3 != p2\r
923 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
924 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
925 Double_t px12 =(photon1+photon2).Px();\r
926 Double_t py12 =(photon1+photon2).Py();\r
927 Double_t pz12 =(photon1+photon2).Pz();\r
928 Double_t e12 =(photon1+photon2).E();\r
929 TLorentzVector photon12(px12, py12,pz12,e12);\r
930 TLorentzVector photon3(p3->Px(),p3->Py(),p3->Pz(),p3->E());\r
931 GetAsyDistPtM(photon12, photon3,fAsy,fDist,fPt,fMass);\r
932\r
933 for(Int_t ipid=0;ipid<fNPID;ipid++){ //pid\r
934 if(p1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
935 p2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
936 p3->IsPIDOK(ipid,AliCaloPID::kPhoton)){\r
937 fMixBPi0GammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module, with pid\r
938 }\r
939 for(Int_t jmod=0;jmod<fNMod;jmod++){ //mod\r
940 for(Int_t kasy=0;kasy<fNAsy;kasy++){ //asy\r
941 if(IsPhotonSelected(p1,ipid,fModCut[jmod]) &&\r
942 IsPhotonSelected(p2,ipid,fModCut[jmod]) &&\r
943 IsPhotonSelected(p3,ipid,fModCut[jmod])){ //photon selectin\r
944 //fill the hist\r
945 for(Int_t mhist=0;mhist<fNHistos;mhist++){ //hist\r
946 Double_t pt1 =mhist *fPerPtBin;\r
947 Double_t pt2 =(mhist+1)*fPerPtBin;\r
948 Int_t index = ipid*fNMod*fNAsy*fNHistos+ jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
949 if(fMass<fInvMassCut && f2gammaAsy<fAsyCut[kasy] && \r
950 fPt<pt2 && fPt>pt1&& fAsy<fAsyCut[kasy] &&fAsy&&fDist&&fPt&&fMass)\r
951 fMixOmegaB[index]->Fill(fMass);\r
952 }//hist \r
953 } //photon selection\r
954 }//asy\r
955 }//mod\r
956 }//pid\r
957 } \r
958 }//ptc3 \r
959 } //pi0 candidate\r
960 } //mix B \r
961\r
962 }//ie\r
963\r
964}\r
965\r
966//______________________________________________________________________________ \r
967void AliAnaNeutralMeson::MixOmegaC(TClonesArray * array1, TClonesArray * array2, TClonesArray * array3)\r
968{\r
969 //omega background\r
970 //three omega background, we classify it A, B and C\r
971 // --A\r
972 // (r1_event1+r2_event1)+r3_event2\r
973 //\r
974 // --B\r
975 // (r1_event1+r2_event2)+r3_event2\r
976 // \r
977 // --C\r
978 // (r1_event1+r2_event2)+r3_event3\r
979 Int_t nclusters1 = array1->GetEntries();\r
980 Int_t nclusters2 = array2->GetEntries();\r
981 Int_t nclusters3 = array3->GetEntries();\r
982 for(Int_t ie=0;ie<nclusters1;ie++){\r
983 for(Int_t je=0;je<nclusters2;je++){\r
984 AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array1->At(ie);\r
985 AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array2->At(je);\r
986 Double_t f2gammaAsy;\r
987 Get2GammaAsyDistPtM(p1,p2,f2gammaAsy,fDist,fPt,fMass);\r
988 if(TMath::Abs(fMass-fPi0Mass)<fPi0MassPeakWidthCut){\r
989 for(Int_t ke=0;ke<nclusters3;ke++){ \r
990 AliAODPWG4Particle * p3 = (AliAODPWG4Particle *)array3->At(ke); //third photon\r
991 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
992 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
993 Double_t px12 =(photon1+photon2).Px();\r
994 Double_t py12 =(photon1+photon2).Py();\r
995 Double_t pz12 =(photon1+photon2).Pz();\r
996 Double_t e12 =(photon1+photon2).E();\r
997 TLorentzVector photon12(px12, py12,pz12,e12);\r
998 TLorentzVector photon3(p3->Px(),p3->Py(),p3->Pz(),p3->E());\r
999 GetAsyDistPtM(photon12, photon3,fAsy,fDist,fPt,fMass);\r
1000\r
1001 for(Int_t ipid=0;ipid<fNPID;ipid++){\r
1002 if(p1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
1003 p2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
1004 p3->IsPIDOK(ipid,AliCaloPID::kPhoton)){\r
1005 fMixCPi0GammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module, with pid\r
1006 }\r
1007 for(Int_t jmod=0;jmod<fNMod;jmod++){\r
1008 for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
1009 if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && \r
1010 IsPhotonSelected(p2,ipid,fModCut[jmod])&&\r
1011 IsPhotonSelected(p3,ipid,fModCut[jmod])){ \r
1012 \r
1013 //filling the histo\r
1014 for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
1015 Double_t pt1 =mhist *fPerPtBin;\r
1016 Double_t pt2 =(mhist+1)*fPerPtBin;\r
1017 Int_t index = ipid*fNMod*fNAsy*fNHistos+ jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
1018 if(fMass<fInvMassCut && f2gammaAsy<fAsyCut[kasy] && \r
1019 fPt<pt2 && fPt>pt1&& fAsy<fAsyCut[kasy] &&fAsy&&fDist&&fPt&&fMass)\r
1020 fMixOmegaC[index]->Fill(fMass);\r
1021 }//hist\r
1022 }//photon selection\r
1023 }//asy\r
1024 }//mod\r
1025 }//pid\r
1026 }//ke\r
1027 }//pi0 candidate\r
1028 } //je\r
1029 }//ie\r
1030\r
1031} \r
1032 \r
1033//______________________________________________________________________________\r
1034void AliAnaNeutralMeson::PhotonAcceptance(AliStack * stack)\r
1035{\r
1036 //photon Acceptance\r
1037 for(Int_t i=0 ; i<stack->GetNtrack(); i++){\r
1038 TParticle * prim = stack->Particle(i) ;\r
1039 Int_t pdg = prim->GetPdgCode() ;\r
1040// Int_t ndau = prim->GetNDaughters();\r
1041// Int_t iphot1=prim->GetFirstDaughter() ;\r
1042// Int_t iphot2=prim->GetLastDaughter() ;\r
1043 if(pdg ==22 && stack->IsPhysicalPrimary(i)) {\r
1044 Double_t photonPt = prim->Pt() ;\r
1045 //printf("photon, pt %2.2f\n",photonPt); \r
1046 Double_t photonY = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ; \r
1047 Double_t phi = TMath::RadToDeg()*prim->Phi() ;\r
1048 if(TMath::Abs(photonY) < 0.5){ fhPrimPhotonPt->Fill(photonPt); }\r
1049 fhPrimPhotonY ->Fill(photonY) ;\r
1050 fhPrimPhotonPhi->Fill(phi) ;\r
1051 //Check if both photons hit Calorimeter\r
1052 Bool_t inacceptance = kFALSE;\r
1053#ifdef __PHOSGEO__\r
1054 Int_t mod1 ;\r
1055 Double_t x1,z1;\r
1056 if(fCalorimeter == "PHOS" && fPHOSGeo->ImpactOnEmc(prim,mod1,z1,x1) ) inacceptance = kTRUE;\r
1057 //printf("In REAL PHOS acceptance? %d\n",inacceptance);\r
1058#else\r
1059 TLorentzVector lv1;\r
1060 prim->Momentum(lv1);\r
1061 if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) ) inacceptance = kTRUE ;\r
1062 //printf("In %s fFidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);\r
1063#endif\r
1064 if(inacceptance){\r
1065 fhPrimPhotonAccPt->Fill(photonPt) ;\r
1066 fhPrimPhotonAccPhi->Fill(phi) ;\r
1067 fhPrimPhotonAccY->Fill(photonY) ;\r
1068 }//Accepted\r
1069 \r
1070 }// Primary photon\r
1071 }//loop on primaries \r
1072}\r
1073 \r
1074//______________________________________________________________________________\r
1075void AliAnaNeutralMeson::Pi0EtaAcceptance(AliStack * stack)\r
1076{\r
1077 //pi0 and eta Acceptance\r
1078 for(Int_t i=0 ; i<stack->GetNprimary(); i++){\r
1079 TParticle * prim = stack->Particle(i) ;\r
1080 Int_t pdg = prim->GetPdgCode() ;\r
1081 if(pdg==111 || pdg==221){\r
1082 Double_t pi0Pt = prim->Pt() ;\r
1083 //printf(" pi0, pt %2.2f\n",pi0Pt);\r
1084 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;\r
1085 Double_t phi = TMath::RadToDeg()*prim->Phi();\r
1086 if(pdg==111){\r
1087 if(TMath::Abs(pi0Y) < 0.5) fhPrimPi0Pt->Fill(pi0Pt) ;\r
1088 fhPrimPi0Y ->Fill(pi0Y) ;\r
1089 fhPrimPi0Phi->Fill(phi) ;\r
1090 } \r
1091 else{\r
1092 if(TMath::Abs(pi0Y) < 0.5) fhPrimEtaPt->Fill(pi0Pt) ;\r
1093 fhPrimEtaY ->Fill(pi0Y) ;\r
1094 fhPrimEtaPhi->Fill(phi) ;\r
1095 }\r
1096 //Check if both photons hit Calorimeter\r
1097 Int_t ndau = prim->GetNDaughters();\r
1098 Int_t iphot1=prim->GetFirstDaughter() ;\r
1099 Int_t iphot2=prim->GetLastDaughter() ;\r
1100 if(ndau==2 && ndau>1 && iphot1>-1 && iphot1<stack->GetNtrack() \r
1101 && iphot2>-1 && iphot2<stack->GetNtrack()){ \r
1102 TParticle * phot1 = stack->Particle(iphot1) ;\r
1103 TParticle * phot2 = stack->Particle(iphot2) ;\r
1104 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){\r
1105 Bool_t inacceptance = kFALSE; \r
1106#ifdef __PHOSGEO__\r
1107 Int_t mod, mod1,mod2 ; \r
1108 Double_t x, z, x1,z1,x2,z2;\r
1109 if(fCalorimeter == "PHOS" && fPHOSGeo->ImpactOnEmc(prim, mod, z,x) && fPHOSGeo->ImpactOnEmc(phot1,mod1,z1,x1) \r
1110 && fPHOSGeo->ImpactOnEmc(phot2,mod2,z2,x2)) inacceptance = kTRUE;\r
1111 \r
1112 //printf("In REAL PHOS acceptance? %d\n",inacceptance);\r
1113#else \r
1114 TLorentzVector lv0, lv1, lv3;\r
1115 prim->Momentum(lv0);\r
1116 phot1->Momentum(lv1);\r
1117 phot2->Momentum(lv3);\r
1118 if(GetFidutialCut()->IsInFidutialCut(lv0,fCalorimeter) \r
1119 &&GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) \r
1120 && GetFidutialCut()->IsInFidutialCut(lv3,fCalorimeter)) inacceptance = kTRUE ;\r
1121#endif\r
1122// printf("In %s fFidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);\r
1123 if(inacceptance){\r
1124 if(pdg==111){\r
1125 fhPrimPi0AccPt->Fill(pi0Pt) ;\r
1126 fhPrimPi0AccPhi->Fill(phi) ; fhPrimPi0AccY->Fill(pi0Y) ;\r
1127 }\r
1128 else{\r
1129 fhPrimEtaAccPt->Fill(pi0Pt) ; fhPrimEtaAccPhi->Fill(phi) ;\r
1130 fhPrimEtaAccY->Fill(pi0Y) ;\r
1131 }\r
1132 }//Accepted\r
1133 }// 2 photons \r
1134 }//Check daughters exist\r
1135 }// Primary pi0 \r
1136 }//loop on primaries \r
1137}\r
1138 \r
1139//______________________________________________________________________________\r
1140void AliAnaNeutralMeson::OmegaAcceptance(AliStack * stack)\r
1141{\r
1142 //acceptance for omega->pi0+gamma \r
1143 for(Int_t i=0 ; i<stack->GetNprimary(); i++){ //aaa\r
1144 TParticle * prim = stack->Particle(i) ;\r
1145 Int_t pdg = prim->GetPdgCode();\r
1146 if(pdg == 223){ //bbb\r
1147 Double_t pt = prim->Pt() ;\r
1148 Double_t y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;\r
1149 Double_t phi = TMath::RadToDeg()*prim->Phi() ;\r
1150 if(TMath::Abs(y) < 0.5){ fhPrimOmegaPt->Fill(pt) ; }\r
1151 fhPrimOmegaY ->Fill(y) ; \r
1152 fhPrimOmegaPhi->Fill(phi) ;\r
1153 //check if omega->pi0+gamma->(gamma+gamma)+gamma hit Calorimeter\r
1154 Int_t ndau = prim->GetNDaughters();\r
1155 Int_t iphot1=prim->GetFirstDaughter() ;\r
1156 Int_t iphot2=prim->GetLastDaughter() ;\r
1157 if(ndau==2 && iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){ //ccc\r
1158 TParticle * phot1 = stack->Particle(iphot1) ;\r
1159 TParticle * phot2 = stack->Particle(iphot2) ;\r
1160 Int_t pdg1 = phot1->GetPdgCode();\r
1161 Int_t pdg2 = phot2->GetPdgCode(); \r
1162 Int_t ndau2 = phot2->GetNDaughters();\r
1163 if(ndau2==2 && phot1 && pdg1==22 &&phot2 && pdg2==111){ //ddd\r
1164 Int_t jphot1 = phot2->GetFirstDaughter() ;\r
1165 Int_t jphot2 = phot2->GetLastDaughter() ;\r
1166 if(jphot1>-1 && jphot1<stack->GetNtrack() &&\r
1167 jphot2>-1 && jphot2<stack->GetNtrack()){ //eee \r
1168 TParticle * phot3 = stack->Particle(jphot1) ;\r
1169 TParticle * phot4 = stack->Particle(jphot2) ;\r
1170 Int_t pdg3 = phot3->GetPdgCode();\r
1171 Int_t pdg4 = phot4->GetPdgCode();\r
1172 if(phot3 && pdg3==22 && phot4 && pdg4==22){ //fff\r
1173 Bool_t inacceptance = kFALSE;\r
1174#ifdef __PHOSGEO__\r
1175 Int_t mod, mod1,mod2,mod3 ;\r
1176 Double_t x,z, x1,z1,x2,z2, x3,z3 ;\r
1177 //3 gammas hit Calorimeter\r
1178 if(fCalorimeter == "PHOS" && fPHOSGeo->ImpactOnEmc(prim,mod,z,x) \r
1179 && fPHOSGeo->ImpactOnEmc(phot1,mod1,z1,x1) \r
1180 && fPHOSGeo->ImpactOnEmc(phot3,mod2,z2,x2) \r
1181 && fPHOSGeo->ImpactOnEmc(phot4,mod3,z3,x3))\r
1182 inacceptance = kTRUE;\r
1183 \r
1184#else\r
1185 TLorentzVector lv0,lv1, lv2, lv3; \r
1186 prim->Momentum(lv0);\r
1187 phot1->Momentum(lv1);\r
1188 phot3->Momentum(lv2);\r
1189 phot4->Momentum(lv3);\r
1190 if(GetFidutialCut()->IsInFidutialCut(lv0,fCalorimeter) &&\r
1191 GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && \r
1192 GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter) && \r
1193 GetFidutialCut()->IsInFidutialCut(lv3,fCalorimeter)) \r
1194 inacceptance = kTRUE ;\r
1195#endif\r
1196 if(inacceptance){ //ggg\r
1197 fhPrimOmegaAccPt->Fill(pt) ;\r
1198 fhPrimOmegaAccPhi->Fill(phi) ;\r
1199 fhPrimOmegaAccY->Fill(y) ;\r
1200 } //ggg\r
1201 }//fff\r
1202 }//eee\r
1203 }//ddd\r
1204 }//ccc\r
1205 } //bbb\r
1206 }//aaa \r
1207} \r
1208 \r
1209void AliAnaNeutralMeson::GetMCDistAsy(TParticle *p1,TParticle *p2,Double_t &dist,Double_t &asy)\r
1210{ \r
1211 //get dist and asy for MC particle \r
1212 Double_t eta1 = p1->Eta(); \r
1213 Double_t eta2 = p2->Eta(); \r
1214 Double_t phi1 = p1->Phi();\r
1215 Double_t phi2 = p2->Phi(); \r
1216 Double_t e1 = p1->Energy(); \r
1217 Double_t e2 = p2->Energy(); \r
1218 dist = TMath::Sqrt( (eta1-eta2)*(eta1-eta2) + (phi1-phi2)*(phi1-phi2));\r
1219 if(e1+e2) asy = TMath::Abs(e1-e2)/(e1+e2);\r
1220} \r
1221\r
1222//______________________________________________________________________________\r
1223//void AliAnaNeutralMeson::Terminate(TList * outList) \r
1224//{\r
1225// //Do some calculations and plots from the final histograms.\r
1226// printf("AliAnaNeutralMeson::Terminate() - Not implemented yet\n");\r
1227//\r
1228//}\r
1229\r
1230\r