]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGUD/UPC/AliAnalysisTaskUpcPsi2s.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bugfix AddTask macro
[u/mrichter/AliRoot.git] / PWGUD / UPC / AliAnalysisTaskUpcPsi2s.cxx
CommitLineData
3d16cd00 1/*************************************************************************
2* Copyright(c) 1998-2008, 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// c++ headers
17#include <iostream>
18#include <string.h>
19
20// root headers
21#include "TH1I.h"
22#include "TTree.h"
23#include "TClonesArray.h"
24#include "TParticle.h"
25#include "TObjString.h"
26#include "TFile.h"
17c65770 27#include "TDatabasePDG.h"
28#include "TLorentzVector.h"
3d16cd00 29
30// aliroot headers
31#include "AliAnalysisManager.h"
32#include "AliInputEventHandler.h"
33#include "AliESDEvent.h"
34#include "AliAODEvent.h"
35#include "AliMCEvent.h"
36#include "AliAODVZERO.h"
37#include "AliAODZDC.h"
38#include "AliESDVZERO.h"
39#include "AliESDZDC.h"
40#include "AliPIDResponse.h"
41#include "AliAODTrack.h"
42#include "AliAODPid.h"
43#include "AliAODVertex.h"
44#include "AliESDVertex.h"
45#include "AliMultiplicity.h"
46#include "AliESDtrack.h"
47#include "AliESDMuonTrack.h"
48#include "AliAODMCParticle.h"
49#include "AliMCParticle.h"
e63dc644 50#include "AliCentrality.h"
489951cd 51#include "AliKFVertex.h"
52#include "AliExternalTrackParam.h"
8a2486ba 53#include "AliTriggerAnalysis.h"
3d16cd00 54
55// my headers
56#include "AliAnalysisTaskUpcPsi2s.h"
57
58ClassImp(AliAnalysisTaskUpcPsi2s);
59
b8f92f9d 60using std::cout;
61using std::endl;
62
3d16cd00 63//trees for UPC analysis,
64// michal.broz@cern.ch
65
66//_____________________________________________________________________________
67AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s()
8a2486ba 68 : AliAnalysisTaskSE(),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
d457ec9e 69 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
03cff2a0 70 fTOFtrig1(0), fTOFtrig2(0),
489951cd 71 fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
72 fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
3d16cd00 73 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
74 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
8a2486ba 75 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
4e649ec8 76 fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),
77 fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
f052ef6f 78 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
79 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
17c65770 80
3d16cd00 81{
82
83//Dummy constructor
84
85}//AliAnalysisTaskUpcPsi2s
86
87
88//_____________________________________________________________________________
89AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
8a2486ba 90 : AliAnalysisTaskSE(name),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
d457ec9e 91 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
03cff2a0 92 fTOFtrig1(0), fTOFtrig2(0),
489951cd 93 fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
94 fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
3d16cd00 95 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
96 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
8a2486ba 97 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
4e649ec8 98 fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),
99 fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
f052ef6f 100 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
101 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
17c65770 102
3d16cd00 103{
104
105 // Constructor
106 if( strstr(name,"ESD") ) fType = 0;
107 if( strstr(name,"AOD") ) fType = 1;
108
109 Init();
110
111 DefineOutput(1, TTree::Class());
112 DefineOutput(2, TTree::Class());
f052ef6f 113 DefineOutput(3, TList::Class());
17c65770 114 DefineOutput(4, TList::Class());
3d16cd00 115
116}//AliAnalysisTaskUpcPsi2s
117
118//_____________________________________________________________________________
119void AliAnalysisTaskUpcPsi2s::Init()
120{
121
122 for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
8a2486ba 123 for(Int_t i=0; i<4; i++) {
124 fTOFphi[i] = -666;
125 fPIDMuon[i] = -666;
126 fPIDElectron[i] = -666;
127 fPIDPion[i] = -666;
128 fTriggerInputsMC[i] = kFALSE;
129 }
489951cd 130 for(Int_t i=0; i<3; i++){
131 fVtxPos[i] = -666;
132 fVtxErr[i] = -666;
133 fKfVtxPos[i] = -666;
134 }
3d16cd00 135
136}//Init
137
138//_____________________________________________________________________________
139AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
140{
141 // Destructor
142 if(fJPsiTree){
143 delete fJPsiTree;
144 fJPsiTree = 0x0;
145 }
146 if(fPsi2sTree){
147 delete fPsi2sTree;
148 fPsi2sTree = 0x0;
149 }
f052ef6f 150 if(fListTrig){
151 delete fListTrig;
152 fListTrig = 0x0;
3d16cd00 153 }
46e1d1dc 154 if(fListHist){
155 delete fListHist;
156 fListHist = 0x0;
17c65770 157 }
3d16cd00 158
159}//~AliAnalysisTaskUpcPsi2s
160
161
162//_____________________________________________________________________________
163void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
164{
8a2486ba 165 //PID response
166 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
167 AliInputEventHandler *inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
168 fPIDResponse = inputHandler->GetPIDResponse();
169
3d16cd00 170 //input file
171 fDataFilnam = new TObjString();
172 fDataFilnam->SetString("");
173
174 //tracks
175 fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
176 fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
177 fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
178 fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
8a2486ba 179 fGenPart = new TClonesArray("TParticle", 1000);
3d16cd00 180
181 //output tree with JPsi candidate events
182 fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
183 fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
184 fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
185 fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
186
187 fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
188 fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
189 fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
190 fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
191 fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
489951cd 192 fJPsiTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
3d16cd00 193 fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
d457ec9e 194
03cff2a0
MB		 195 fJPsiTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
196 fJPsiTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
8b17ae4a 197 fJPsiTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
03cff2a0 198
8a2486ba 199 fJPsiTree ->Branch("fPIDMuon", &fPIDMuon[0], "fPIDMuon[4]/D");
200 fJPsiTree ->Branch("fPIDElectron", &fPIDElectron[0], "fPIDElectron[4]/D");
201 fJPsiTree ->Branch("fPIDPion", &fPIDPion[0], "fPIDPion[4]/D");
202
489951cd 203 fJPsiTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
204 fJPsiTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
d457ec9e
M		 205 fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
206 fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
207
489951cd 208 fJPsiTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
209
3d16cd00 210 fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
211 fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
212 fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
213 fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
214 fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
215 fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
216 fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
217 if( fType == 0 ) {
218 fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
219 }
220 if( fType == 1 ) {
221 fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
222 }
8a2486ba 223 if(isMC) {
224 fJPsiTree ->Branch("fGenPart", &fGenPart);
225 fJPsiTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
226 }
227
3d16cd00 228
229 //output tree with Psi2s candidate events
230 fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
231 fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
232 fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
233 fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
234
235 fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
236 fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
237 fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
238 fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
239 fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
489951cd 240 fPsi2sTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
3d16cd00 241 fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
d457ec9e 242
03cff2a0
MB		 243 fPsi2sTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
244 fPsi2sTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
8b17ae4a 245 fPsi2sTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
03cff2a0 246
8a2486ba 247 fPsi2sTree ->Branch("fPIDMuon", &fPIDMuon[0], "fPIDMuon[4]/D");
248 fPsi2sTree ->Branch("fPIDElectron", &fPIDElectron[0], "fPIDElectron[4]/D");
249 fPsi2sTree ->Branch("fPIDPion", &fPIDPion[0], "fPIDPion[4]/D");
250
489951cd 251 fPsi2sTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
252 fPsi2sTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
d457ec9e
M		 253 fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
254 fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
255
489951cd 256 fPsi2sTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
257
3d16cd00 258 fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
259 fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
260 fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
261 fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
262 fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
263 fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
264 fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
265 if( fType == 0 ) {
266 fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
267 }
268 if( fType == 1 ) {
269 fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
270 }
8a2486ba 271 if(isMC) {
272 fPsi2sTree ->Branch("fGenPart", &fGenPart);
273 fPsi2sTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
274 }
275
3d16cd00 276
f052ef6f 277 fListTrig = new TList();
278 fListTrig ->SetOwner();
279
4e649ec8 280 fHistCcup4TriggersPerRun = new TH1D("fHistCcup4TriggersPerRun", "fHistCcup4TriggersPerRun", 33000, 167000.5, 200000.5);
281 fListTrig->Add(fHistCcup4TriggersPerRun);
f052ef6f 282
4e649ec8 283 fHistCcup7TriggersPerRun = new TH1D("fHistCcup7TriggersPerRun", "fHistCcup7TriggersPerRun", 33000, 167000.5, 200000.5);
284 fListTrig->Add(fHistCcup7TriggersPerRun);
285
286 fHistCcup2TriggersPerRun = new TH1D("fHistCcup2TriggersPerRun", "fHistCcup2TriggersPerRun", 33000, 167000.5, 200000.5);
287 fListTrig->Add(fHistCcup2TriggersPerRun);
288
289 fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 33000, 167000.5, 200000.5);
f052ef6f 290 fListTrig->Add(fHistZedTriggersPerRun);
291
4e649ec8 292 fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 33000, 167000.5, 200000.5);
f052ef6f 293 fListTrig->Add(fHistCvlnTriggersPerRun);
294
4e649ec8 295 fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 33000, 167000.5, 200000.5);
e63dc644
MB		 296 fListTrig->Add(fHistMBTriggersPerRun);
297
4e649ec8 298 fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 33000, 167000.5, 200000.5);
e63dc644
MB		 299 fListTrig->Add(fHistCentralTriggersPerRun);
300
4e649ec8 301 fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 33000, 167000.5, 200000.5);
e63dc644
MB		 302 fListTrig->Add(fHistSemiCentralTriggersPerRun);
303
46e1d1dc 304 fListHist = new TList();
305 fListHist ->SetOwner();
306
d1fa9007 307 TString CutNameJPsi[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Two good tracks",
17c65770 308 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
d1fa9007 309 fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",13,0.5,13.5);
310 for (Int_t i = 0; i<13; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
46e1d1dc 311 fListHist->Add(fHistNeventsJPsi);
17c65770 312
313 fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
46e1d1dc 314 fListHist->Add(fHistTPCsignalJPsi);
17c65770 315
316 fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
46e1d1dc 317 fListHist->Add(fHistDiLeptonPtJPsi);
17c65770 318
319 fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
320 fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
46e1d1dc 321 fListHist->Add(fHistDiElectronMass);
17c65770 322
323 fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
324 fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
46e1d1dc 325 fListHist->Add(fHistDiMuonMass);
17c65770 326
d1fa9007 327 TString CutNamePsi2s[14] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Four good tracks",
17c65770 328 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
329
d1fa9007 330 fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",14,0.5,14.5);
331 for (Int_t i = 0; i<14; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
46e1d1dc 332 fListHist->Add(fHistNeventsPsi2s);
17c65770 333
334 fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
335 fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
336 fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
46e1d1dc 337 fListHist->Add(fHistPsi2sMassVsPt);
17c65770 338
339 fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
340 fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
46e1d1dc 341 fListHist->Add(fHistPsi2sMassCoherent);
17c65770 342
3d16cd00 343 PostData(1, fJPsiTree);
344 PostData(2, fPsi2sTree);
f052ef6f 345 PostData(3, fListTrig);
46e1d1dc 346 PostData(4, fListHist);
3d16cd00 347
348}//UserCreateOutputObjects
349
350//_____________________________________________________________________________
351void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
352{
353
354 //cout<<"#################### Next event ##################"<<endl;
355
f052ef6f 356 if( fType == 0 ){
357 RunESDtrig();
358 if(fRunHist) RunESDhist();
359 if(fRunTree) RunESDtree();
360 }
361
362 if( fType == 1 ){
363 RunAODtrig();
17c65770 364 if(fRunHist) RunAODhist();
365 if(fRunTree) RunAODtree();
366 }
3d16cd00 367
368}//UserExec
17c65770 369//_____________________________________________________________________________
f052ef6f 370void AliAnalysisTaskUpcPsi2s::RunAODtrig()
371{
372
373 //input event
374 AliAODEvent *aod = (AliAODEvent*) InputEvent();
375 if(!aod) return;
376
377 fRunNum = aod ->GetRunNumber();
378 //Trigger
379 TString trigger = aod->GetFiredTriggerClasses();
380
4e649ec8 381 if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
382 if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
383 if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
f052ef6f 384
e63dc644
MB		 385 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
386 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
f052ef6f 387
388 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
389 if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
e63dc644
MB		 390
391 //MB, Central and SemiCentral triggers
392 AliCentrality *centrality = aod->GetCentrality();
393 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
394
03cff2a0
MB		 395 Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
396 //Double_t percentile = centrality->GetCentralityPercentile("V0M");
e63dc644 397
d1fa9007 398 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
399
4e649ec8 400 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
d1fa9007 401
402 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
3c79e118 403
f052ef6f 404PostData(3, fListTrig);
405
406}
407//_____________________________________________________________________________
17c65770 408void AliAnalysisTaskUpcPsi2s::RunAODhist()
409{
410
411 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
412
413 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
414 Double_t muonMass = partMuon->Mass();
415
416 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
417 Double_t electronMass = partElectron->Mass();
418
419 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
420 Double_t pionMass = partPion->Mass();
421
422 //input event
423 AliAODEvent *aod = (AliAODEvent*) InputEvent();
424 if(!aod) return;
425
426 fHistNeventsJPsi->Fill(1);
427 fHistNeventsPsi2s->Fill(1);
428
429 //Trigger
430 TString trigger = aod->GetFiredTriggerClasses();
431
8a2486ba 432 if(!isMC && !trigger.Contains("CCUP") ) return;
17c65770 433
434 fHistNeventsJPsi->Fill(2);
435 fHistNeventsPsi2s->Fill(2);
436
437 //primary vertex
438 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
439 fVtxContrib = fAODVertex->GetNContributors();
440 if(fVtxContrib < 2) return;
441
442 fHistNeventsJPsi->Fill(3);
443 fHistNeventsPsi2s->Fill(3);
444
17c65770 445 //VZERO, ZDC
446 AliAODVZERO *fV0data = aod ->GetVZEROData();
f052ef6f 447 AliAODZDC *fZDCdata = aod->GetZDCData();
17c65770 448
449 fV0Adecision = fV0data->GetV0ADecision();
450 fV0Cdecision = fV0data->GetV0CDecision();
451 if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
452
d1fa9007 453 fHistNeventsJPsi->Fill(4);
454 fHistNeventsPsi2s->Fill(4);
455
f052ef6f 456 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
457 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
458
e63dc644 459 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
f052ef6f 460
d1fa9007 461 fHistNeventsJPsi->Fill(5);
462 fHistNeventsPsi2s->Fill(5);
17c65770 463
17c65770 464 //Two tracks loop
fc8ccb59 465 Int_t nGoodTracks = 0;
17c65770 466 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
467
468 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
469 Short_t qLepton[4], qPion[4];
470 UInt_t nLepton=0, nPion=0, nHighPt=0;
489951cd 471 Double_t fRecTPCsignal[5];
17c65770 472 Int_t mass[3]={-1,-1,-1};
73c6b208 473
fc8ccb59 474
73c6b208 475 //Four track loop
476 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
477 AliAODTrack *trk = aod->GetTrack(itr);
478 if( !trk ) continue;
d1fa9007 479 if(!(trk->TestFilterBit(1<<0))) continue;
73c6b208 480
481 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
482 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
483 if(trk->GetTPCNcls() < 50)continue;
484 if(trk->Chi2perNDF() > 4)continue;
485 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
fc8ccb59
MB		 486 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
487 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
488 delete trk_clone;
489
73c6b208 490 if(TMath::Abs(dca[1]) > 2) continue;
491
492 TrackIndex[nGoodTracks] = itr;
493 nGoodTracks++;
494
495 if(nGoodTracks > 4) break;
496 }//Track loop
497
17c65770 498 nLepton=0; nPion=0; nHighPt=0;
499 mass[0]= -1; mass[1]= -1, mass[2]= -1;
500
501 if(nGoodTracks == 4){
d1fa9007 502 fHistNeventsPsi2s->Fill(6);
17c65770 503 for(Int_t i=0; i<4; i++){
504 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
505
506 if(trk->Pt() > 1){
489951cd 507 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
17c65770 508 qLepton[nLepton] = trk->Charge();
489951cd 509 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
17c65770 510 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
511 mass[nLepton] = 0;
512 }
489951cd 513 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
17c65770 514 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
515 mass[nLepton] = 1;
516 }
517 nLepton++;
518 }
519 else{
520 qPion[nPion] = trk->Charge();
521 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
522 nPion++;
523 }
524
525 if(nLepton > 2 || nPion > 2) break;
526 }
527 if((nLepton == 2) && (nPion == 2)){
d1fa9007 528 fHistNeventsPsi2s->Fill(7);
529 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(8);
530 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(9);
531 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
17c65770 532 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
d1fa9007 533 fHistNeventsPsi2s->Fill(11);
17c65770 534 if(mass[0] == mass[1]) {
d1fa9007 535 fHistNeventsPsi2s->Fill(12);
17c65770 536 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
537 vDilepton = vLepton[0]+vLepton[1];
538 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
539 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
d1fa9007 540 if(mass[0] == 0) fHistNeventsPsi2s->Fill(13);
541 if(mass[0] == 1) fHistNeventsPsi2s->Fill(14);
17c65770 542 }
543 }
544 }
545 }
546
fc8ccb59
MB		 547 nGoodTracks = 0;
548 //Two track loop
549 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
550 AliAODTrack *trk = aod->GetTrack(itr);
551 if( !trk ) continue;
d1fa9007 552 if(!(trk->TestFilterBit(1<<0))) continue;
fc8ccb59
MB		 553
554 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
555 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
556 if(trk->GetTPCNcls() < 70)continue;
557 if(trk->Chi2perNDF() > 4)continue;
558 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
559 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
560 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
561 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
562 delete trk_clone;
563 if(TMath::Abs(dca[1]) > 2) continue;
564 if(TMath::Abs(dca[0]) > 0.2) continue;
565
566 TrackIndex[nGoodTracks] = itr;
567 nGoodTracks++;
568
569 if(nGoodTracks > 2) break;
570 }//Track loop
571
572 nLepton=0; nPion=0; nHighPt=0;
573 mass[0]= -1; mass[1]= -1, mass[2]= -1;
574
575 if(nGoodTracks == 2){
d1fa9007 576 fHistNeventsJPsi->Fill(6);
fc8ccb59
MB		 577 for(Int_t i=0; i<2; i++){
578 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
579 if(trk->Pt() > 1) nHighPt++;
489951cd 580 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
fc8ccb59 581 qLepton[nLepton] = trk->Charge();
489951cd 582 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
fc8ccb59
MB		 583 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
584 mass[nLepton] = 0;
585 }
489951cd 586 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
fc8ccb59
MB		 587 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
588 mass[nLepton] = 1;
589 }
590 nLepton++;
591 }
592 if(nLepton == 2){
d1fa9007 593 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(7);
fc8ccb59 594 if(qLepton[0]*qLepton[1] < 0){
d1fa9007 595 fHistNeventsJPsi->Fill(8);
fc8ccb59 596 if(nHighPt > 0){
d1fa9007 597 fHistNeventsJPsi->Fill(9);
489951cd 598 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
d1fa9007 599 if(nHighPt == 2) fHistNeventsJPsi->Fill(10);
fc8ccb59 600 if(mass[0] == mass[1] && mass[0] != -1) {
d1fa9007 601 fHistNeventsJPsi->Fill(11);
fc8ccb59
MB		 602 vCandidate = vLepton[0]+vLepton[1];
603 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
604 if(mass[0] == 0) {
605 fHistDiMuonMass->Fill(vCandidate.M());
d1fa9007 606 fHistNeventsJPsi->Fill(12);
fc8ccb59
MB		 607 }
608 if(mass[0] == 1) {
609 fHistDiElectronMass->Fill(vCandidate.M());
d1fa9007 610 fHistNeventsJPsi->Fill(13);
fc8ccb59
MB		 611 }
612 }
613 }
614 }
615 }
616 }
617
618
46e1d1dc 619 PostData(4, fListHist);
17c65770 620
621}
3d16cd00 622
623//_____________________________________________________________________________
17c65770 624void AliAnalysisTaskUpcPsi2s::RunAODtree()
3d16cd00 625{
626 //input event
627 AliAODEvent *aod = (AliAODEvent*) InputEvent();
628 if(!aod) return;
629
8a2486ba 630 if(isMC) RunAODMC(aod);
631
3d16cd00 632 //input data
633 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
634 fDataFilnam->Clear();
635 fDataFilnam->SetString(filnam);
636 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
637 fRunNum = aod ->GetRunNumber();
638
3d16cd00 639 //Trigger
640 TString trigger = aod->GetFiredTriggerClasses();
641
8b17ae4a
MB		 642 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
643 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
644 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
645
646 Bool_t isTriggered = kFALSE;
647 for(Int_t i=0; i<ntrg; i++) {
648 if( fTrigger[i] ) isTriggered = kTRUE;
649 }
8a2486ba 650 if(!isMC && !isTriggered ) return;
3d16cd00 651
652 //trigger inputs
653 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
8b17ae4a 654 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
3d16cd00 655
656 //Event identification
657 fPerNum = aod ->GetPeriodNumber();
658 fOrbNum = aod ->GetOrbitNumber();
659 fBCrossNum = aod ->GetBunchCrossNumber();
660
661 //primary vertex
662 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
663 fVtxContrib = fAODVertex->GetNContributors();
489951cd 664 fVtxPos[0] = fAODVertex->GetX();
665 fVtxPos[1] = fAODVertex->GetY();
666 fVtxPos[2] = fAODVertex->GetZ();
d457ec9e
M		 667 Double_t CovMatx[6];
668 fAODVertex->GetCovarianceMatrix(CovMatx);
489951cd 669 fVtxErr[0] = CovMatx[0];
670 fVtxErr[1] = CovMatx[1];
671 fVtxErr[2] = CovMatx[2];
d457ec9e
M		 672 fVtxChi2 = fAODVertex->GetChi2();
673 fVtxNDF = fAODVertex->GetNDF();
3d16cd00 674
675 //Tracklets
676 fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
677
678 //VZERO, ZDC
679 AliAODVZERO *fV0data = aod ->GetVZEROData();
680 AliAODZDC *fZDCdata = aod->GetZDCData();
681
682 fV0Adecision = fV0data->GetV0ADecision();
683 fV0Cdecision = fV0data->GetV0CDecision();
684 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
685 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
3d16cd00 686
489951cd 687 fNLooseTracks = 0;
688
689 //Track loop - loose cuts
690 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
691 AliAODTrack *trk = aod->GetTrack(itr);
692 if( !trk ) continue;
d1fa9007 693 if(!(trk->TestFilterBit(1<<0))) continue;
489951cd 694
695 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
696 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
697 if(trk->GetTPCNcls() < 20)continue;
698 fNLooseTracks++;
699 }//Track loop -loose cuts
700
17c65770 701 Int_t nGoodTracks=0;
3d16cd00 702 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
fc8ccb59 703
489951cd 704 //Two track loop
3d16cd00 705 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
706 AliAODTrack *trk = aod->GetTrack(itr);
707 if( !trk ) continue;
d1fa9007 708 if(!(trk->TestFilterBit(1<<0))) continue;
709
489951cd 710 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
711 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
8b17ae4a 712 if(trk->GetTPCNcls() < 70)continue;
3d16cd00 713 if(trk->Chi2perNDF() > 4)continue;
8b17ae4a 714 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
3d16cd00 715 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
fc8ccb59
MB		 716 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
717 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
718 delete trk_clone;
3d16cd00 719 if(TMath::Abs(dca[1]) > 2) continue;
38677cd7 720 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
721 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
722
8b17ae4a 723
3d16cd00 724 TrackIndex[nGoodTracks] = itr;
725 nGoodTracks++;
726
8b17ae4a 727 if(nGoodTracks > 2) break;
3d16cd00 728 }//Track loop
8a2486ba 729
730 fJPsiAODTracks->Clear("C");
8b17ae4a 731 if(nGoodTracks == 2){
489951cd 732
733 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
734 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
735 Double_t muonMass = partMuon->Mass();
736 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
737 Double_t electronMass = partElectron->Mass();
738 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
739 Double_t pionMass = partPion->Mass();
740
741 Double_t KFcov[21];
742 Double_t KFpar[6];
743 Double_t KFmass = pionMass;
744 Double_t fRecTPCsignal;
745 AliKFParticle *KFpart[2];
746 AliKFVertex *KFvtx = new AliKFVertex();
747 KFvtx->SetField(aod->GetMagneticField());
748
8b17ae4a 749 for(Int_t i=0; i<2; i++){
3d16cd00 750 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
751
752 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
fc8ccb59
MB		 753 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
754 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
755 delete trk_clone;
8b17ae4a
MB		 756
757 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
758 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
3d16cd00 759
8a2486ba 760 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
761 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
762 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
763
489951cd 764 trk->GetPosition(KFpar);
765 trk->PxPyPz(KFpar+3);
766 trk->GetCovMatrix(KFcov);
767
768 if(trk->Pt() > 1){
769 fRecTPCsignal = trk->GetTPCsignal();
770 if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
771 if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
772 }
773 else KFmass = pionMass;
774
775 KFpart[i] = new AliKFParticle();
776 KFpart[i]->SetField(aod->GetMagneticField());
777 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
778 KFvtx->AddDaughter(*KFpart[i]);
779
780
8b17ae4a 781 Double_t pos[3]={0,0,0};
489951cd 782 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
783 parTrk->CopyFromVTrack((AliVTrack*) trk);
784 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
8b17ae4a
MB		 785 else {
786 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
787 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
489951cd 788 }
789 delete parTrk;
3d16cd00 790 }
489951cd 791 fKfVtxPos[0]= KFvtx->GetX();
792 fKfVtxPos[1]= KFvtx->GetY();
793 fKfVtxPos[2]= KFvtx->GetZ();
794 for(UInt_t i=0; i<2; i++)delete KFpart[i];
795 delete KFvtx;
796
8a2486ba 797 if(!isMC) fJPsiTree ->Fill();
3d16cd00 798 }
8b17ae4a
MB		 799
800 nGoodTracks = 0;
801 //Four track loop
73c6b208 802 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
803 AliAODTrack *trk = aod->GetTrack(itr);
804 if( !trk ) continue;
d1fa9007 805 if(!(trk->TestFilterBit(1<<0))) continue;
73c6b208 806
489951cd 807 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
808 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
8b17ae4a 809 if(trk->GetTPCNcls() < 50)continue;
73c6b208 810 if(trk->Chi2perNDF() > 4)continue;
811 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
fc8ccb59
MB		 812 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
813 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
814 delete trk_clone;
8b17ae4a 815 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
73c6b208 816 if(TMath::Abs(dca[1]) > 2) continue;
38677cd7 817 Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
818 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
8b17ae4a 819
73c6b208 820 TrackIndex[nGoodTracks] = itr;
821 nGoodTracks++;
822
8b17ae4a 823 if(nGoodTracks > 4) break;
73c6b208 824 }//Track loop
489951cd 825
8a2486ba 826 fPsi2sAODTracks->Clear("C");
8b17ae4a 827 if(nGoodTracks == 4){
d1fa9007 828
489951cd 829 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
830 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
831 Double_t muonMass = partMuon->Mass();
832 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
833 Double_t electronMass = partElectron->Mass();
834 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
835 Double_t pionMass = partPion->Mass();
836
837 Double_t KFcov[21];
838 Double_t KFpar[6];
839 Double_t KFmass = pionMass;
840 Double_t fRecTPCsignal;
841 AliKFParticle *KFpart[4];
842 AliKFVertex *KFvtx = new AliKFVertex();
843 KFvtx->SetField(aod->GetMagneticField());
844
8b17ae4a 845 for(Int_t i=0; i<4; i++){
3d16cd00 846 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
847
848 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
fc8ccb59
MB		 849 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
850 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
851 delete trk_clone;
3d16cd00 852
8b17ae4a
MB		 853 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
854 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
8a2486ba 855
856
857 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
858 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
859 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
8b17ae4a 860
489951cd 861 trk->GetPosition(KFpar);
862 trk->PxPyPz(KFpar+3);
863 trk->GetCovMatrix(KFcov);
864
865 if(trk->Pt() > 1){
866 fRecTPCsignal = trk->GetTPCsignal();
867 if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
868 if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
869 }
870 else KFmass = pionMass;
871
872 KFpart[i] = new AliKFParticle();
873 KFpart[i]->SetField(aod->GetMagneticField());
874 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
875 KFvtx->AddDaughter(*KFpart[i]);
876
8b17ae4a 877 Double_t pos[3]={0,0,0};
489951cd 878 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
879 parTrk->CopyFromVTrack((AliVTrack*) trk);
880 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
8b17ae4a
MB		 881 else {
882 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
883 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
489951cd 884 }
885 delete parTrk;
3d16cd00 886 }
489951cd 887 fKfVtxPos[0]= KFvtx->GetX();
888 fKfVtxPos[1]= KFvtx->GetY();
889 fKfVtxPos[2]= KFvtx->GetZ();
890 for(UInt_t i=0; i<4; i++)delete KFpart[i];
891 delete KFvtx;
8a2486ba 892 if(!isMC) fPsi2sTree ->Fill();
893 }
894
895 if(isMC){
896 fJPsiTree ->Fill();
897 fPsi2sTree ->Fill();
3d16cd00 898 }
03cff2a0
MB		 899
900 PostData(1, fJPsiTree);
901 PostData(2, fPsi2sTree);
3d16cd00 902
903}//RunAOD
904
8a2486ba 905
906//_____________________________________________________________________________
907void AliAnalysisTaskUpcPsi2s::RunAODMC(AliAODEvent *aod)
908{
909
910 fGenPart->Clear("C");
911
912 TClonesArray *arrayMC = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
913 if(!arrayMC) return;
914
915 Int_t nmc=0;
916 //loop over mc particles
917 for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
918 AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(imc);
919 if(!mcPart) continue;
920
921 if(mcPart->GetMother() >= 0) continue;
922
923 TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
924 part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
925 part->SetPdgCode(mcPart->GetPdgCode());
926 part->SetUniqueID(imc);
927 }//loop over mc particles
928
929}//RunAODMC
930
931
3d16cd00 932//_____________________________________________________________________________
f052ef6f 933void AliAnalysisTaskUpcPsi2s::RunESDtrig()
934{
935
936 //input event
937 AliESDEvent *esd = (AliESDEvent*) InputEvent();
938 if(!esd) return;
939
940 fRunNum = esd ->GetRunNumber();
941 //Trigger
942 TString trigger = esd->GetFiredTriggerClasses();
943
4e649ec8 944 if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
945 if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
946 if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
f052ef6f 947
e63dc644
MB		 948 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
949 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
f052ef6f 950
951 if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
952
03cff2a0
MB		 953 //MB, Central and SemiCentral triggers
954 AliCentrality *centrality = esd->GetCentrality();
955 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
956
957 //Double_t percentile = centrality->GetCentralityPercentile("V0M");
958 Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
959
d1fa9007 960 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
03cff2a0 961
4e649ec8 962 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
03cff2a0 963
d1fa9007 964 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
03cff2a0
MB		 965
966
f052ef6f 967PostData(3, fListTrig);
968
969}
970//_____________________________________________________________________________
971void AliAnalysisTaskUpcPsi2s::RunESDhist()
972{
973
974 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
975
976 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
977 Double_t muonMass = partMuon->Mass();
978
979 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
980 Double_t electronMass = partElectron->Mass();
981
982 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
983 Double_t pionMass = partPion->Mass();
984
985 //input event
986 AliESDEvent *esd = (AliESDEvent*) InputEvent();
987 if(!esd) return;
988
989 fHistNeventsJPsi->Fill(1);
990 fHistNeventsPsi2s->Fill(1);
991
992 //Trigger
993 TString trigger = esd->GetFiredTriggerClasses();
994
8a2486ba 995 if(!isMC && !trigger.Contains("CCUP") ) return;
f052ef6f 996
997 fHistNeventsJPsi->Fill(2);
998 fHistNeventsPsi2s->Fill(2);
999
1000 //primary vertex
1001 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1002 fVtxContrib = fESDVertex->GetNContributors();
1003 if(fVtxContrib < 2) return;
1004
1005 fHistNeventsJPsi->Fill(3);
1006 fHistNeventsPsi2s->Fill(3);
1007
1008 //VZERO, ZDC
1009 AliESDVZERO *fV0data = esd->GetVZEROData();
1010 AliESDZDC *fZDCdata = esd->GetESDZDC();
1011
1012 fV0Adecision = fV0data->GetV0ADecision();
1013 fV0Cdecision = fV0data->GetV0CDecision();
1014 if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
1015
1016 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1017 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
03cff2a0 1018 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
f052ef6f 1019
1020 fHistNeventsJPsi->Fill(4);
1021 fHistNeventsPsi2s->Fill(4);
1022
1023 Int_t nGoodTracks=0;
1024 //Two tracks loop
1025 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1026
1027 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1028 Short_t qLepton[4], qPion[4];
1029 UInt_t nLepton=0, nPion=0, nHighPt=0;
489951cd 1030 Double_t fRecTPCsignal[5];
f052ef6f 1031 Int_t mass[3]={-1,-1,-1};
1032
03cff2a0 1033 //Two Track loop
f052ef6f 1034 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1035 AliESDtrack *trk = esd->GetTrack(itr);
1036 if( !trk ) continue;
1037
1038 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1039 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
03cff2a0 1040 if(trk->GetTPCNcls() < 70)continue;
f052ef6f 1041 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
03cff2a0 1042 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
f052ef6f 1043 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1044 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1045 trk->GetImpactParameters(dca[0],dca[1]);
1046 if(TMath::Abs(dca[1]) > 2) continue;
03cff2a0 1047 if(TMath::Abs(dca[1]) > 0.2) continue;
f052ef6f 1048
1049 TrackIndex[nGoodTracks] = itr;
1050 nGoodTracks++;
03cff2a0 1051 if(nGoodTracks > 2) break;
f052ef6f 1052 }//Track loop
03cff2a0 1053
f052ef6f 1054
1055 if(nGoodTracks == 2){
1056 fHistNeventsJPsi->Fill(5);
1057 for(Int_t i=0; i<2; i++){
1058 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1059 if(trk->Pt() > 1) nHighPt++;
489951cd 1060 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
f052ef6f 1061 qLepton[nLepton] = trk->Charge();
489951cd 1062 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
f052ef6f 1063 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1064 mass[nLepton] = 0;
1065 }
489951cd 1066 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
f052ef6f 1067 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1068 mass[nLepton] = 1;
1069 }
1070 nLepton++;
1071 }
1072 if(nLepton == 2){
1073 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
1074 if(qLepton[0]*qLepton[1] < 0){
1075 fHistNeventsJPsi->Fill(7);
1076 if(nHighPt > 0){
1077 fHistNeventsJPsi->Fill(8);
489951cd 1078 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
f052ef6f 1079 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
1080 if(mass[0] == mass[1] && mass[0] != -1) {
1081 fHistNeventsJPsi->Fill(10);
1082 vCandidate = vLepton[0]+vLepton[1];
1083 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
1084 if(mass[0] == 0) {
1085 fHistDiMuonMass->Fill(vCandidate.M());
1086 fHistNeventsJPsi->Fill(11);
1087 }
1088 if(mass[0] == 1) {
1089 fHistDiElectronMass->Fill(vCandidate.M());
1090 fHistNeventsJPsi->Fill(12);
1091 }
1092 }
1093 }
1094 }
1095 }
1096 }
03cff2a0 1097 nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
f052ef6f 1098 mass[0]= -1; mass[1]= -1, mass[2]= -1;
1099
03cff2a0
MB		 1100 //Four Track loop
1101 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1102 AliESDtrack *trk = esd->GetTrack(itr);
1103 if( !trk ) continue;
1104
1105 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1106 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1107 if(trk->GetTPCNcls() < 50)continue;
1108 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1109 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1110 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1111 trk->GetImpactParameters(dca[0],dca[1]);
1112 if(TMath::Abs(dca[1]) > 2) continue;
1113
1114 TrackIndex[nGoodTracks] = itr;
1115 nGoodTracks++;
1116 if(nGoodTracks > 4) break;
1117 }//Track loop
1118
f052ef6f 1119 if(nGoodTracks == 4){
1120 fHistNeventsPsi2s->Fill(5);
1121 for(Int_t i=0; i<4; i++){
1122 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1123
1124 if(trk->Pt() > 1){
489951cd 1125 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
f052ef6f 1126 qLepton[nLepton] = trk->Charge();
489951cd 1127 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
f052ef6f 1128 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1129 mass[nLepton] = 0;
1130 }
489951cd 1131 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
f052ef6f 1132 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1133 mass[nLepton] = 1;
1134 }
1135 nLepton++;
1136 }
1137 else{
1138 qPion[nPion] = trk->Charge();
1139 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1140 nPion++;
1141 }
1142
1143 if(nLepton > 2 || nPion > 2) break;
1144 }
1145 if((nLepton == 2) && (nPion == 2)){
1146 fHistNeventsPsi2s->Fill(6);
1147 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1148 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1149 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1150 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1151 fHistNeventsPsi2s->Fill(10);
1152 if(mass[0] == mass[1]) {
1153 fHistNeventsPsi2s->Fill(11);
1154 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1155 vDilepton = vLepton[0]+vLepton[1];
1156 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1157 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1158 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1159 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1160 }
1161 }
1162 }
1163 }
1164
1165 PostData(4, fListHist);
1166
1167}
1168
1169//_____________________________________________________________________________
1170void AliAnalysisTaskUpcPsi2s::RunESDtree()
3d16cd00 1171{
1172
1173 //input event
1174 AliESDEvent *esd = (AliESDEvent*) InputEvent();
1175 if(!esd) return;
8a2486ba 1176
1177 if(isMC) RunESDMC(esd);
3d16cd00 1178
1179 //input data
1180 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1181 fDataFilnam->Clear();
1182 fDataFilnam->SetString(filnam);
1183 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1184 fRunNum = esd->GetRunNumber();
1185
8b17ae4a 1186 //Trigger
3d16cd00 1187 TString trigger = esd->GetFiredTriggerClasses();
1188
8b17ae4a
MB		 1189 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1190 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
1191 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1192
1193 Bool_t isTriggered = kFALSE;
1194 for(Int_t i=0; i<ntrg; i++) {
1195 if( fTrigger[i] ) isTriggered = kTRUE;
1196 }
8a2486ba 1197 if(!isMC && !isTriggered ) return;
8b17ae4a 1198
3d16cd00 1199 //trigger inputs
1200 fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1201 fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
03cff2a0
MB		 1202
1203 //TOF trigger info (0OMU)
1204 fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1205 fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
3d16cd00 1206
1207 //Event identification
1208 fPerNum = esd->GetPeriodNumber();
1209 fOrbNum = esd->GetOrbitNumber();
1210 fBCrossNum = esd->GetBunchCrossNumber();
1211
1212 //primary vertex
1213 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1214 fVtxContrib = fESDVertex->GetNContributors();
489951cd 1215 fVtxPos[0] = fESDVertex->GetX();
1216 fVtxPos[1] = fESDVertex->GetY();
1217 fVtxPos[2] = fESDVertex->GetZ();
03cff2a0
MB		 1218 Double_t CovMatx[6];
1219 fESDVertex->GetCovarianceMatrix(CovMatx);
489951cd 1220 fVtxErr[0] = CovMatx[0];
1221 fVtxErr[1] = CovMatx[1];
1222 fVtxErr[2] = CovMatx[2];
03cff2a0
MB		 1223 fVtxChi2 = fESDVertex->GetChi2();
1224 fVtxNDF = fESDVertex->GetNDF();
3d16cd00 1225
1226 //Tracklets
1227 fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1228
1229 //VZERO, ZDC
1230 AliESDVZERO *fV0data = esd->GetVZEROData();
1231 AliESDZDC *fZDCdata = esd->GetESDZDC();
1232
1233 fV0Adecision = fV0data->GetV0ADecision();
1234 fV0Cdecision = fV0data->GetV0CDecision();
1235 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1236 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
489951cd 1237
1238 fNLooseTracks = 0;
1239
1240 //Track loop - loose cuts
1241 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1242 AliESDtrack *trk = esd->GetTrack(itr);
1243 if( !trk ) continue;
1244
1245 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1246 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1247 if(trk->GetTPCNcls() < 20)continue;
1248 fNLooseTracks++;
1249 }//Track loop -loose cuts
3d16cd00 1250
17c65770 1251 Int_t nGoodTracks=0;
3d16cd00 1252 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
17c65770 1253
03cff2a0 1254 //Two Track loop
3d16cd00 1255 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1256 AliESDtrack *trk = esd->GetTrack(itr);
1257 if( !trk ) continue;
1258
1259 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1260 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
03cff2a0 1261 if(trk->GetTPCNcls() < 70)continue;
3d16cd00 1262 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
03cff2a0 1263 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
3d16cd00 1264 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1265 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1266 trk->GetImpactParameters(dca[0],dca[1]);
1267 if(TMath::Abs(dca[1]) > 2) continue;
03cff2a0 1268 if(TMath::Abs(dca[1]) > 0.2) continue;
3d16cd00 1269
1270 TrackIndex[nGoodTracks] = itr;
1271 nGoodTracks++;
03cff2a0 1272 if(nGoodTracks > 2) break;
3d16cd00 1273 }//Track loop
1274
8a2486ba 1275 fJPsiESDTracks->Clear("C");
3d16cd00 1276 if(nGoodTracks == 2){
1277 for(Int_t i=0; i<2; i++){
1278 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1279
1280 AliExternalTrackParam cParam;
1281 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1282
1283 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
1284
8a2486ba 1285 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1286 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1287 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1288
8b17ae4a
MB		 1289 Double_t pos[3]={0,0,0};
1290 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1291 else {
1292 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1293 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1294 }
3d16cd00 1295 }
8a2486ba 1296 if(!isMC) fJPsiTree ->Fill();
3d16cd00 1297 }
1298
03cff2a0
MB		 1299 nGoodTracks = 0;
1300 //Four track loop
1301 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1302 AliESDtrack *trk = esd->GetTrack(itr);
1303 if( !trk ) continue;
1304
1305 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1306 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1307 if(trk->GetTPCNcls() < 50)continue;
1308 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1309 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1310 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1311 trk->GetImpactParameters(dca[0],dca[1]);
1312 if(TMath::Abs(dca[1]) > 2) continue;
1313
1314 TrackIndex[nGoodTracks] = itr;
1315 nGoodTracks++;
1316 if(nGoodTracks > 4) break;
1317 }//Track loop
1318
8a2486ba 1319 fPsi2sESDTracks->Clear("C");
3d16cd00 1320 if(nGoodTracks == 4){
1321 for(Int_t i=0; i<4; i++){
1322 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1323
1324 AliExternalTrackParam cParam;
1325 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1326
8b17ae4a 1327 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
8a2486ba 1328
1329 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1330 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1331 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
8b17ae4a
MB		 1332
1333 Double_t pos[3]={0,0,0};
1334 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1335 else {
1336 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1337 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1338 }
3d16cd00 1339 }
8a2486ba 1340 if(!isMC) fPsi2sTree ->Fill();
3d16cd00 1341 }
03cff2a0 1342
8a2486ba 1343 if(isMC){
1344 fJPsiTree ->Fill();
1345 fPsi2sTree ->Fill();
1346 }
1347
03cff2a0
MB		 1348 PostData(1, fJPsiTree);
1349 PostData(2, fPsi2sTree);
3d16cd00 1350
1351}//RunESD
1352
8a2486ba 1353
1354//_____________________________________________________________________________
1355void AliAnalysisTaskUpcPsi2s::RunESDMC(AliESDEvent* esd)
1356{
1357
1358 AliTriggerAnalysis *fTrigAna = new AliTriggerAnalysis();
1359 fTrigAna->SetAnalyzeMC(isMC);
1360
1361 if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) > 1) fTriggerInputsMC[0] = kTRUE;
1362 if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) < 2) fTriggerInputsMC[0] = kFALSE;
1363
1364 fTriggerInputsMC[1] = esd->GetHeader()->IsTriggerInputFired("0OMU");
1365 fTriggerInputsMC[2] = esd->GetHeader()->IsTriggerInputFired("0VBA");
1366 fTriggerInputsMC[3] = esd->GetHeader()->IsTriggerInputFired("0VBC");
1367
1368 fGenPart->Clear("C");
1369
1370 AliMCEvent *mc = MCEvent();
1371 if(!mc) return;
1372
1373 Int_t nmc = 0;
1374 //loop over mc particles
1375 for(Int_t imc=0; imc<mc->GetNumberOfTracks(); imc++) {
1376 AliMCParticle *mcPart = (AliMCParticle*) mc->GetTrack(imc);
1377 if(!mcPart) continue;
1378
1379 if(mcPart->GetMother() >= 0) continue;
1380
1381 TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1382 part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1383 part->SetPdgCode(mcPart->PdgCode());
1384 part->SetUniqueID(imc);
1385 }//loop over mc particles
1386
1387}//RunESDMC
1388
1389
1390
3d16cd00 1391//_____________________________________________________________________________
1392void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
1393{
1394
1395 cout<<"Analysis complete."<<endl;
1396}//Terminate
1397