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