]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/PHOSTasks/PHOS_PbPb/AliAnalysisTaskPi0FlowMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changed Input::Input() to work with AnalysisResult.root files with pPb data.
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / AliAnalysisTaskPi0FlowMC.cxx
CommitLineData
444647ad 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// Extension to Pi0FLOw, mimicing AliPHOSHijingEfficiency
17// by Dmitri Peressounko, 05.02.2013
18// Authors: Henrik Qvigstad, Dmitri Peressounko
19// Date : 05.04.2013
20/* $Id$ */
21
22#include "TChain.h"
23#include "TTree.h"
24#include "TObjArray.h"
25#include "TF1.h"
26#include "TFile.h"
27#include "TH1F.h"
28#include "TH2F.h"
29#include "TH2I.h"
30#include "TH3F.h"
31#include "TParticle.h"
32#include "TCanvas.h"
33#include "TStyle.h"
34#include "TRandom.h"
35#include "TROOT.h"
36#include "THashList.h"
37
38
39#include "AliAnalysisManager.h"
40#include "AliMCEventHandler.h"
41#include "AliMCEvent.h"
42#include "AliStack.h"
43#include "AliAnalysisTaskSE.h"
44#include "AliPHOSHijingEfficiency.h"
45#include "AliCaloPhoton.h"
46#include "AliPHOSGeometry.h"
47#include "AliPHOSEsdCluster.h"
48#include "AliPHOSCalibData.h"
49#include "AliESDEvent.h"
50#include "AliESDCaloCells.h"
51#include "AliESDVertex.h"
52#include "AliESDtrackCuts.h"
53#include "AliLog.h"
54#include "AliPID.h"
55#include "AliCDBManager.h"
56#include "AliCentrality.h"
57#include "AliESDtrackCuts.h"
58#include "AliEventplane.h"
59#include "TProfile.h"
60#include <TPDGCode.h>
61#include "AliOADBContainer.h"
62
63
64#include "AliAnalysisTaskPi0FlowMC.h"
65
66ClassImp(AliAnalysisTaskPi0FlowMC);
67
68//TODO: rnlin?
69
70AliAnalysisTaskPi0FlowMC::AliAnalysisTaskPi0FlowMC(const char* name, AliAnalysisTaskPi0Flow::Period period)
71: AliAnalysisTaskPi0Flow(name, period),
72 fStack(0x0)
73{
74}
75
76AliAnalysisTaskPi0FlowMC::~AliAnalysisTaskPi0FlowMC()
77{
78}
79
80void AliAnalysisTaskPi0FlowMC::MakeMCHistograms()
81{
82 //AliAnalysisTaskPi0Flow::MakeMCHistograms();
83
84 // MC Generated histograms
85 char key[55];
86 for(Int_t cent=0; cent < fCentEdges.GetSize()-1; cent++){
87 snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
88 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
89 snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
90 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
91 snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
92 fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ;
93 snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
94 fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
95 snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
96 fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
97 snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
98 fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
99 snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
100 fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
101 snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
102 fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
103 snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
104 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
105 snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
106 fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
107 snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
108 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
109 snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
110 fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
111 }
112 fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
113 fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
114 fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
115
116
117 Int_t nPt = 200;
118 Double_t ptMin = 0;
119 Double_t ptMax = 20;
120 fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.));
121 fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.));
122 fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.));
123 fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.));
124 fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.));
125 fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.));
126// fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.));
127// fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.));
128
129 const Int_t nM = 500;
130 const Double_t mMin = 0.0;
131 const Double_t mMax = 1.0;
132
133 for(Int_t cen=0; cen < fCentEdges.GetSize()-1; cen++){
134 fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
135 fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
136 fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
137 fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
138 fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
139 fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
140 fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
141 fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
142 fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
143 fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
144 fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
145 fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
146 fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
147
148 //pairs from common parents
149 fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
150 fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
151 fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
152 fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
153 fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
154 fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
155 fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
156
157 //common parent - pi0
158 fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
159 fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
160 fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
161 fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
162 fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
163 fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
164 fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
165 fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
166 fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
167
168 }
169
170
171 //Photon contaminations
172 fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
173 fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
174 fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
175 fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
176 fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.));
177
178 const Int_t nTypes=24 ;
179 char partTypes[nTypes][55] ;
180 snprintf(partTypes[0],55,"hGammaNoPrim") ; //
181 snprintf(partTypes[1],55,"hGammaPhot") ; //
182 snprintf(partTypes[2],55,"hGammaEl") ; //
183 snprintf(partTypes[3],55,"hGammaPi0") ; //
184 snprintf(partTypes[4],55,"hGammaEta") ; //
185 snprintf(partTypes[5],55,"hhGammaOmega") ; //
186 snprintf(partTypes[6],55,"hGammaPipm") ; //
187 snprintf(partTypes[7],55,"hGammaP") ; //
188 snprintf(partTypes[8],55,"hGammaPbar") ; //
189 snprintf(partTypes[9],55,"hGammaN") ; //
190 snprintf(partTypes[10],55,"hGammaNbar") ; //
191 snprintf(partTypes[11],55,"hGammaK0S") ; //
192 snprintf(partTypes[12],55,"hGammaK0L") ; //
193 snprintf(partTypes[13],55,"hGammaKpm") ; //
194 snprintf(partTypes[14],55,"hGammaKstar") ; //
195 snprintf(partTypes[15],55,"hGammaDelta") ; //
196 snprintf(partTypes[16],55,"hGammaOtherCharged") ; //
197 snprintf(partTypes[17],55,"hGammaOtherNeutral") ; //
198 snprintf(partTypes[18],55,"hGammaPipmGamma") ; //
199 snprintf(partTypes[19],55,"hGammaPipmEl") ; //
200 snprintf(partTypes[20],55,"hGammaPipmOther") ; //
201 snprintf(partTypes[21],55,"hGammaPipmDirect") ; //
202 snprintf(partTypes[22],55,"hGammaPipmp") ; //
203 snprintf(partTypes[23],55,"hGammaPipmn") ; //
204
205 const Int_t nPID=12 ;
545b989c 206 char cPID[12][25] ;
444647ad 207 snprintf(cPID[0],25,"All") ;
208 snprintf(cPID[1],25,"Allcore") ;
209 snprintf(cPID[2],25,"CPV") ;
210 snprintf(cPID[3],25,"CPVcore") ;
211 snprintf(cPID[4],25,"CPV2") ;
212 snprintf(cPID[5],25,"CPV2core") ;
213 snprintf(cPID[6],25,"Disp") ;
214 snprintf(cPID[7],25,"Dispcore") ;
215 snprintf(cPID[8],25,"Disp2") ;
216 snprintf(cPID[9],25,"Disp2core") ;
217 snprintf(cPID[10],25,"Both") ;
218 snprintf(cPID[11],25,"Bothcore") ;
219
220 for(Int_t itype=0; itype<nTypes; itype++){
221 for(Int_t iPID=0; iPID<nPID; iPID++){
222 for(Int_t cen=0; cen<5; cen++){
223 fOutputContainer->Add(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax));
224 }
225 }
226 }
227}
228
229void AliAnalysisTaskPi0FlowMC::DoMC()
230{
231 //AliAnalysisTaskPi0Flow::DoMC();
232 fStack = GetMCStack();
233
234 //TODO: Geometry IHEP?
235 //TODO: decalib.?
236 //TODO: PHOS matrix?
237 //TODO: Centrality.
238
239 FillMCHist();
240
241 // Proccess all selected clusters
242 for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
243 //Bool_t sure= kTRUE;
244 //Int_t primary=FindPrimary(clu,sure) ; //номер праймари частицы в стеке
245 //ph->SetPrimary(primary) ;
246 //ph->SetWeight(PrimaryWeight(primary)) ;
247 }
248
249 FillSecondaries() ;
250}
251
252
253
254//___________________________________________________________________________
255void AliAnalysisTaskPi0FlowMC::FillMCHist(){
256 //fill histograms for efficiensy etc. calculation
257
258 //---------First pi0/eta-----------------------------
259 char partName[10] ;
260 char hkey[55] ;
261
262 if(!fStack) return ;
263 for(Int_t i=0;i<fStack->GetNtrack();i++){
264 TParticle* particle = fStack->Particle(i);
265 if(particle->GetPdgCode() == kPi0)
266 snprintf(partName,10,"pi0") ;
267 else
268 if(particle->GetPdgCode() == kEta)
269 snprintf(partName,10,"eta") ;
270 else
271 if(particle->GetPdgCode() == kGamma)
272 snprintf(partName,10,"gamma") ;
273 else
274 continue ;
275
276 //Primary particle
277 Double_t r=particle->R() ;
278 Double_t pt = particle->Pt() ;
279 //Distribution over vertex
280 FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
281
282 if(r >kRCut)
283 continue ;
284
285 //Total number of pi0 with creation radius <1 cm
286 Double_t weight = PrimaryParticleWeight(particle) ;
287 snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCentBin) ;
288 FillHistogram(hkey,pt,weight) ;
289 if(TMath::Abs(particle->Y())<0.12){
290 snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCentBin) ;
291 FillHistogram(hkey,pt,weight) ;
292 }
293
294 snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCentBin) ;
295 FillHistogram(hkey,particle->Y(),weight) ;
296
297 Double_t phi=particle->Phi() ;
298 while(phi<0.)phi+=TMath::TwoPi() ;
299 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
300 snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCentBin) ;
301 FillHistogram(hkey,phi,weight) ;
302 }
303}
304
305//________________________________________________________________________
306void AliAnalysisTaskPi0FlowMC::FillSecondaries(){
307 //Sort secondaires
308
309 //Fill spectra of primary particles
310 //with proper weight
311 std::cout << fStack << std::endl;
312 AliInfo("start");
313 for(Int_t i=0; i<fStack->GetNtrack(); i++){
314 TParticle * p = fStack->Particle(i) ;
315 if(p->R()>kRCut)
316 continue ;
317 if(TMath::Abs(p->Y())>0.5)
318 continue ;
319 Double_t w = PrimaryParticleWeight(p) ;
320 Int_t primPdgCode=p->GetPdgCode() ;
321 switch(primPdgCode){
322 case kGamma: FillHistogram(Form("hPrimPhot_cen%d",fCentBin),p->Pt(),w);
323 break ;
324 case kElectron:
325 case -kElectron:
326 FillHistogram(Form("hPrimEl_cen%d",fCentBin),p->Pt(),w);
327 break ;
328 case kPi0:
329 FillHistogram(Form("hPrimPi0_cen%d",fCentBin),p->Pt(),w);
330 break ;
331 case kEta:
332 FillHistogram(Form("hPrimEta_cen%d",fCentBin),p->Pt(),w);
333 break ;
334 case kPiPlus:
335 case kPiMinus:
336 FillHistogram(Form("hPrimPipm_cen%d",fCentBin),p->Pt(),w);
337 break ;
338 case kProton: //p
339 FillHistogram(Form("hPrimP_cen%d",fCentBin),p->Pt(),w);
340 break ;
341 case kProtonBar: //pbar
342 FillHistogram(Form("hPrimPbar_cen%d",fCentBin),p->Pt(),w);
343 break ;
344 case kNeutron: //n
345 FillHistogram(Form("hPrimN_cen%d",fCentBin),p->Pt(),w);
346 break ;
347 case kNeutronBar: //nbar
348 FillHistogram(Form("hPrimNbar_cen%d",fCentBin),p->Pt(),w);
349 break ;
350 case 310: //nbar
351 FillHistogram(Form("hPrimK0S_cen%d",fCentBin),p->Pt(),w);
352 break ;
353 case 130: //nbar
354 FillHistogram(Form("hPrimK0L_cen%d",fCentBin),p->Pt(),w);
355 break ;
356 case 321: //K+
357 case -321: //K-
358 FillHistogram(Form("hPrimKpm_cen%d",fCentBin),p->Pt(),w);
359 break ;
360 default: //other
361 FillHistogram(Form("hPrimOther_cen%d",fCentBin),p->Pt(),w);
362 }
363 }
364 AliInfo("Origins of secondary pi0s");
365 //Origins of secondary pi0s
366 for(Int_t i=0; i<fStack->GetNtrack(); i++){
367 TParticle * p = fStack->Particle(i) ;
368 if(p->GetPdgCode()!=111)
369 continue ;
370 FillHistogram("Vertex",p->Pt(),p->R());
371 if(p->R()<kRCut)
372 continue ;
373 Double_t phi=p->Phi() ;
374 while(phi<0.)phi+=TMath::TwoPi() ;
375 while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
376 FillHistogram("hSecondPi0RphiZ",p->R(),phi,p->Vz()) ;
377 Double_t w = PrimaryParticleWeight(p) ;
378 FillHistogram("hSecondPi0RE",p->R(),p->Pt(),w) ;
379 }
380
381 TLorentzVector p1;
382
383 Int_t inPHOS=fCaloPhotonsPHOS->GetEntries() ;
384 for (Int_t i1=0; i1<inPHOS-1; i1++) {
385 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
386 Double_t w1=ph1->GetWeight() ;
387 for (Int_t i2=i1+1; i2<inPHOS; i2++) {
388 AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
389 TLorentzVector p12 = *ph1 + *ph2;
390 Double_t w2=ph2->GetWeight() ;
391 Double_t w = TMath::Sqrt(w1*w2) ;
392 FillHistogram(Form("hParentAll_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
393 Int_t prim=FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary()) ;
394 if(prim>-1){
395 TParticle * particle = (TParticle *)fStack->Particle(prim);
396 FillHistogram("hMass_R",p12.M(),p12.Pt(),TMath::Sqrt(particle->R()*particle->R()+particle->Vz()*particle->Vz())) ;
397
398
399 Int_t pdgCode=particle->GetPdgCode() ;
400 if(pdgCode!=111){ //common parent - not pi0
401 if(pdgCode==22)
402 FillHistogram(Form("hParentGamma_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
403 else{
404 if(pdgCode==11 || pdgCode==-11){
405 FillHistogram(Form("hParentEl_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
406 }
407 else{
408 if(InPi0mass(p12.M() ,p12.Pt())){
409 printf("Common parent: %d \n",pdgCode) ;
410 }
411 FillHistogram(Form("hParentOther_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
412 }
413 }//Not photons
414 }//Common parent not pi0
415 else{ //common parent - pi0
416 FillHistogram(Form("hParentPi0_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
417 FillHistogram(Form("Real_pi_R"),p12.M(),p12.Pt(),particle->R(),w) ;
418 FillHistogram(Form("Real_pi_Z"),p12.M(),p12.Pt(),particle->Vz(),w) ;
419 if(particle->R()<kRCut && TMath::Abs(particle->Vz())<fMaxAbsVertexZ){
420 FillHistogram(Form("hParentDirPi0_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
421 continue ;
422 }
423 //Common particle pi0, created off-vertex
424 Int_t primPi0=particle->GetFirstMother();
425 if(primPi0==-1){
426 FillHistogram(Form("hParentPi0NoPrim_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
427 }
428 else{
429 Int_t primPdgCode=fStack->Particle(primPi0)->GetPdgCode() ;
430 switch(primPdgCode){
431 case 221: FillHistogram(Form("hParentPi0Eta_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //eta
432 break ;
433 case 223: FillHistogram(Form("hParentPi0Omega_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //omega
434 break ;
435 case 211: //pi+-
436 case -211: FillHistogram(Form("hParentPi0Pipm_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
437 break ;
438 case 321: //K+-
439 case -321: FillHistogram(Form("hParentPi0Kpm_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
440 break ;
441 case 310: FillHistogram(Form("hParentPi0Ks_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // K0S
442 break ;
443 case 130: FillHistogram(Form("hParentPi0Kl_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // K0L
444 break ;
445 case 2212: //p
446 case 2112: //n
447 FillHistogram(Form("hParentPi0pn_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // pn
448 break ;
449 case -2212: //pbar
450 case -2112: //nbar
451 FillHistogram(Form("hParentPi0antipn_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // pn
452 break ;
453 default: //other
454 FillHistogram(Form("hParentPi0Other_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
455 }//switch
456 }//pi0 with primary
457 }//common parent - pi0
458 }//there is common primary
459 }//seond photon loop
460 }//first photon loop
461
462
463 //Now look at photon contaiminations
464 for (Int_t i1=0; i1<inPHOS-1; i1++) {
465 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
466 Int_t iprim = ph1->GetPrimary() ;
467 if(iprim<0)
468 FillAllHistograms(Form("hGammaNoPrim_cen%d",fCentBin),ph1) ; //
469 else{
470 //Find primary at vertex
471 TParticle * primPHOS = fStack->Particle(iprim) ;
472 Int_t iprimV=primPHOS->GetFirstMother();
473 TParticle * primVtx = primPHOS ;
474 while((iprimV>-1) && primVtx->R()>kRCut){
475 primVtx = fStack->Particle(iprimV) ;
476 iprimV=primVtx->GetFirstMother();
477 }
478
479 //photon
480 Int_t primPdgCode=primVtx->GetPdgCode() ;
481 switch(primPdgCode){
482 case 22: FillAllHistograms("hGammaPhot",ph1);
483 break ;
484 case 11:
485 case -11:
486 FillAllHistograms("hGammaEl",ph1);
487 break ;
488 case 111:
489 FillAllHistograms("hGammaPi0",ph1);
490 break ;
491 case 221:
492 FillAllHistograms("hGammaEta",ph1);
493 break ;
494 case 223: FillAllHistograms("hGammaOmega",ph1) ; //omega
495 break ;
496 case 211:
497 case -211:
498 FillAllHistograms("hGammaPipm",ph1);
499 //Find particle entered PHOS
500 if(primVtx == primPHOS)
501 FillAllHistograms("hGammaPipmDirect",ph1);
502 else{
503 Int_t primPdgPHOS=primPHOS->GetPdgCode() ;
504 if(primPdgPHOS==22){
505 FillAllHistograms("hGammaPipmGamma",ph1);
506 FillHistogram("hPipmGammaConvR",ph1->Pt(),primPHOS->R());
507 FillHistogram("hPipmGammaConvRZ",primPHOS->Vz(),primPHOS->R());
508 break ;
509 }
510 if(TMath::Abs(primPdgPHOS)==11){
511 FillAllHistograms("hGammaPipmEl",ph1);
512 FillHistogram("hPipmElConvR",ph1->Pt(),primPHOS->R());
513 break ;
514 }
515 if(TMath::Abs(primPdgPHOS)==2212){
516 FillAllHistograms("hGammaPipmp",ph1);
517 FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
518 break ;
519 }
520 if(TMath::Abs(primPdgPHOS)==2112){
521 FillAllHistograms("hGammaPipmn",ph1);
522 FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
523 break ;
524 }
525 FillAllHistograms("hGammaPipmOther",ph1);
526 FillHistogram("hPipmOtherConvR",ph1->Pt(),primPHOS->R());
527 }
528 break ;
529 case 2212: //p
530 FillAllHistograms("hGammaP",ph1);
531 break ;
532 case -2212: //pbar
533 FillAllHistograms("hGammaPbar",ph1);
534 break ;
535 case 2112: //n
536 FillAllHistograms("hGammaN",ph1);
537 break ;
538 case -2112: //nbar
539 FillAllHistograms("hGammaNbar",ph1) ; // pn
540 break ;
541 case 310: //nbar
542 FillAllHistograms("hGammaK0S",ph1) ; // pn
543 break ;
544 case 130: //nbar
545 FillAllHistograms("hGammaK0L",ph1) ; // pn
546 break ;
547 case 321: //K+
548 case -321: //K-
549 FillAllHistograms("hGammaKpm",ph1) ; // pn
550 break ;
551 case -323:
552 case 323:
553 case -313:
554 case 313: FillAllHistograms("hGammaKstar",ph1) ; // K*(892)
555 break ;
556
557 case -2224 : //Deltas
558 case 2224 : //Deltas
559 case -2214 : //Deltas
560 case 2214 : //Deltas
561 case -2114 : //Deltas
562 case 2114 : //Deltas
563 case -1114 : //Deltas
564 case 1114 : //Deltas
565 FillAllHistograms("hGammaDelta",ph1) ; // pn
566 break ;
567 default: //other
568 if(primVtx->GetPDG()->Charge())
569 FillAllHistograms("hGammaOtherCharged",ph1) ; //
570 else
571 FillAllHistograms("hGammaOtherNeutral",ph1) ; //
572 }
573 }
574
575 }//single photons
576}
577
578//_____________________________________________________________________________
579void AliAnalysisTaskPi0FlowMC::FillAllHistograms(const char * particleName,AliCaloPhoton * ph)
580{
581 //Fill All PID histograms
582
583 Double_t w=ph->GetWeight() ;
584 Double_t pt = ph->Pt() ;
585 Double_t ptC=ph->GetMomV2()->Pt() ;
586 FillHistogram(Form("%s_All_cen%d",particleName,fCentBin),pt,w) ;
587 FillHistogram(Form("%s_Allcore_cen%d",particleName,fCentBin),ptC,w) ;
588 if(ph->IsCPVOK()){
589 FillHistogram(Form("%s_CPV_cen%d",particleName,fCentBin),pt,w) ;
590 FillHistogram(Form("%s_CPVcore_cen%d",particleName,fCentBin),ptC,w) ;
591 }
592 if(ph->IsCPV2OK()){
593 FillHistogram(Form("%s_CPV2_cen%d",particleName,fCentBin),pt,w) ;
594 FillHistogram(Form("%s_CPV2core_cen%d",particleName,fCentBin),ptC,w) ;
595 }
596 if(ph->IsDispOK()){
597 FillHistogram(Form("%s_Disp_cen%d",particleName,fCentBin),pt,w) ;
598 FillHistogram(Form("%s_Dispcore_cen%d",particleName,fCentBin),ptC,w) ;
599 if(ph->IsDisp2OK()){
600 FillHistogram(Form("%s_Disp2_cen%d",particleName,fCentBin),pt,w) ;
601 FillHistogram(Form("%s_Disp2core_cen%d",particleName,fCentBin),ptC,w) ;
602 }
603 if(ph->IsCPVOK()){
604 FillHistogram(Form("%s_Both_cen%d",particleName,fCentBin),pt,w) ;
605 FillHistogram(Form("%s_Bothcore_cen%d",particleName,fCentBin),ptC,w) ;
606 }
607 }
608}
609
610
611//___________________________________________________________________________
612Double_t AliAnalysisTaskPi0FlowMC::PrimaryWeight(Int_t primary){
613 //Check who is the primary and introduce weight to correct primary spectrum
614
615 if(primary<0 || primary>=fStack->GetNtrack())
616 return 1 ;
617 //trace primaries up to IP
618 TParticle* particle = fStack->Particle(primary);
619 Double_t r=particle->R() ;
620 Int_t mother = particle->GetFirstMother() ;
621 while(mother>-1){
622 if(r<1. && particle->GetPdgCode()==111)
623 break ;
624 particle = fStack->Particle(mother);
625 mother = particle->GetFirstMother() ;
626 r=particle->R() ;
627 }
628
629 return TMath::Max(0.,PrimaryParticleWeight(particle)) ;
630}
631//________________________________________________________________________
632Double_t AliAnalysisTaskPi0FlowMC::PrimaryParticleWeight(TParticle * particle){
545b989c 633
444647ad 634 Int_t pdg = particle->GetPdgCode() ;
635 Int_t type=0 ;
636 if(pdg == 111 || TMath::Abs(pdg)==211){
637 type =1 ;
638 }
639 else{
640 if(TMath::Abs(pdg)<1000){ //Kaon-like
641 type =2 ;
642 }
643 else
644 type = 3; //baryons
645 }
646
647 Double_t pt = particle->Pt() ;
648 if(type==1){
649 if(fCentBin==0) //0-5
650 return (1.662990+1.140890*pt-0.192088*pt*pt)/(1.-0.806630*pt+0.304771*pt*pt)+0.141690*pt ;
651 if(fCentBin==1) //5-10
652 return (1.474351+0.791492*pt-0.066369*pt*pt)/(1.-0.839338*pt+0.317312*pt*pt)+0.093289*pt ;
653 if(fCentBin==2) //10-20
654 return (1.174728+0.959681*pt-0.137695*pt*pt)/(1.-0.788873*pt+0.299538*pt*pt)+0.128759*pt ;
655 if(fCentBin==3) //20-40
656 return (0.927335+0.475349*pt+0.004364*pt*pt)/(1.-0.817966*pt+0.309787*pt*pt)+0.086899*pt ;
657 if(fCentBin==4) //40-60
658 return (0.676878+0.190680*pt+0.077031*pt*pt)/(1.-0.790623*pt+0.305183*pt*pt)+0.064510*pt ;
659 if(fCentBin==5) //60-80
660 return (0.684726-0.606262*pt+0.409963*pt*pt)/(1.-1.080061*pt+0.456933*pt*pt)+0.005151*pt ;
661 }
662 if(type==2){
663 if(fCentBin==0) //0-5
664 return (-0.417131+2.253936*pt-0.337731*pt*pt)/(1.-0.909892*pt+0.316820*pt*pt)+0.157312*pt ;
665 if(fCentBin==1) //5-10
666 return (-0.352275+1.844466*pt-0.248598*pt*pt)/(1.-0.897048*pt+0.316462*pt*pt)+0.132461*pt ;
667 if(fCentBin==2) //10-20
668 return (-0.475481+1.975108*pt-0.336013*pt*pt)/(1.-0.801028*pt+0.276705*pt*pt)+0.188164*pt ;
669 if(fCentBin==3) //20-40
670 return (-0.198954+1.068789*pt-0.103540*pt*pt)/(1.-0.848354*pt+0.299209*pt*pt)+0.112939*pt ;
671 if(fCentBin==4) //40-60
672 return (-0.111052+0.664041*pt-0.019717*pt*pt)/(1.-0.804916*pt+0.300779*pt*pt)+0.095784*pt ;
673 if(fCentBin==5) //0-5
674 return (0.202788-0.439832*pt+0.564585*pt*pt)/(1.-1.254029*pt+0.679444*pt*pt)+0.016235*pt ;
675 }
676 if(type==3){
677 if(fCentBin==0) //0-5
678 return (-1.312732+2.743568*pt-0.375775*pt*pt)/(1.-0.717533*pt+0.164694*pt*pt)+0.164445*pt ;
679 if(fCentBin==1) //5-10
680 return (-1.229425+2.585889*pt-0.330164*pt*pt)/(1.-0.715892*pt+0.167386*pt*pt)+0.133085*pt ;
681 if(fCentBin==2) //10-20
682 return (-1.135677+2.397489*pt-0.320355*pt*pt)/(1.-0.709312*pt+0.164350*pt*pt)+0.146095*pt ;
683 if(fCentBin==3) //20-40
684 return (-0.889993+1.928263*pt-0.220785*pt*pt)/(1.-0.715991*pt+0.174729*pt*pt)+0.095098*pt ;
685 if(fCentBin==4) //40-60
686 return (-0.539237+1.329118*pt-0.115439*pt*pt)/(1.-0.722906*pt+0.186832*pt*pt)+0.059267*pt ;
687 if(fCentBin==5) //60-80
688 return (-0.518126+1.327628*pt-0.130881*pt*pt)/(1.-0.665649*pt+0.184300*pt*pt)+0.081701*pt ;
689 }
690 return 1. ;
691}
692
693//___________________________________________________________________________
694AliStack* AliAnalysisTaskPi0FlowMC::GetMCStack()
695{
696 fStack = 0;
697 AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
698 if(eventHandler){
699 AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
700 if( mcEventHandler)
701 fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
702 }
703 return fStack;
704}