]> 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
Check on contributors of spd vertex
[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()
80b87448 68 : AliAnalysisTaskSE(),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fRunSystematics(kFALSE),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),
a4c57c19 76 fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),fHistCint1TriggersPerRun(0),fHistC0tvxAndCint1TriggersPerRun(0),
4e649ec8 77 fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
80b87448 78 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),fHistDiLeptonMass(0),
79 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
22b569e5 80 fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0),fListPsi2sLoose(0),fListPsi2sTight(0)
17c65770 81
3d16cd00 82{
83
84//Dummy constructor
85
86}//AliAnalysisTaskUpcPsi2s
87
88
89//_____________________________________________________________________________
90AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
80b87448 91 : AliAnalysisTaskSE(name),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fRunSystematics(kFALSE),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
d457ec9e 92 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
03cff2a0 93 fTOFtrig1(0), fTOFtrig2(0),
489951cd 94 fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
95 fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
3d16cd00 96 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
97 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
8a2486ba 98 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
a4c57c19 99 fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),fHistCint1TriggersPerRun(0),fHistC0tvxAndCint1TriggersPerRun(0),
4e649ec8 100 fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
80b87448 101 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),fHistDiLeptonMass(0),
102 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
22b569e5 103 fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0),fListPsi2sLoose(0),fListPsi2sTight(0)
17c65770 104
3d16cd00 105{
106
107 // Constructor
108 if( strstr(name,"ESD") ) fType = 0;
109 if( strstr(name,"AOD") ) fType = 1;
110
111 Init();
112
113 DefineOutput(1, TTree::Class());
114 DefineOutput(2, TTree::Class());
f052ef6f 115 DefineOutput(3, TList::Class());
17c65770 116 DefineOutput(4, TList::Class());
3d16cd00 117
118}//AliAnalysisTaskUpcPsi2s
119
120//_____________________________________________________________________________
121void AliAnalysisTaskUpcPsi2s::Init()
122{
123
124 for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
8a2486ba 125 for(Int_t i=0; i<4; i++) {
126 fTOFphi[i] = -666;
d9117245 127 fPIDTPCMuon[i] = -666;
128 fPIDTPCElectron[i] = -666;
129 fPIDTPCPion[i] = -666;
130 fPIDTPCKaon[i] = -666;
131 fPIDTPCProton[i] = -666;
132
133 fPIDTOFMuon[i] = -666;
134 fPIDTOFElectron[i] = -666;
135 fPIDTOFPion[i] = -666;
136 fPIDTOFKaon[i] = -666;
137 fPIDTOFProton[i] = -666;
138
8a2486ba 139 fTriggerInputsMC[i] = kFALSE;
140 }
489951cd 141 for(Int_t i=0; i<3; i++){
142 fVtxPos[i] = -666;
143 fVtxErr[i] = -666;
144 fKfVtxPos[i] = -666;
d0071d7b 145 fSpdVtxPos[i] = -666;
489951cd 146 }
3d16cd00 147
148}//Init
149
150//_____________________________________________________________________________
151AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
152{
153 // Destructor
154 if(fJPsiTree){
155 delete fJPsiTree;
156 fJPsiTree = 0x0;
157 }
158 if(fPsi2sTree){
159 delete fPsi2sTree;
160 fPsi2sTree = 0x0;
161 }
f052ef6f 162 if(fListTrig){
163 delete fListTrig;
164 fListTrig = 0x0;
3d16cd00 165 }
46e1d1dc 166 if(fListHist){
167 delete fListHist;
168 fListHist = 0x0;
17c65770 169 }
3d16cd00 170
171}//~AliAnalysisTaskUpcPsi2s
172
173
174//_____________________________________________________________________________
175void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
176{
8a2486ba 177 //PID response
178 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
179 AliInputEventHandler *inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
180 fPIDResponse = inputHandler->GetPIDResponse();
181
3d16cd00 182 //input file
183 fDataFilnam = new TObjString();
184 fDataFilnam->SetString("");
185
186 //tracks
187 fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
188 fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
189 fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
190 fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
8a2486ba 191 fGenPart = new TClonesArray("TParticle", 1000);
3d16cd00 192
193 //output tree with JPsi candidate events
194 fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
195 fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
196 fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
197 fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
198
199 fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
200 fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
201 fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
202 fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
203 fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
489951cd 204 fJPsiTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
3d16cd00 205 fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
d457ec9e 206
03cff2a0
MB		 207 fJPsiTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
208 fJPsiTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
8b17ae4a 209 fJPsiTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
03cff2a0 210
d9117245 211 fJPsiTree ->Branch("fPIDTPCMuon", &fPIDTPCMuon[0], "fPIDTPCMuon[4]/D");
212 fJPsiTree ->Branch("fPIDTPCElectron", &fPIDTPCElectron[0], "fPIDTPCElectron[4]/D");
213 fJPsiTree ->Branch("fPIDTPCPion", &fPIDTPCPion[0], "fPIDTPCPion[4]/D");
214 fJPsiTree ->Branch("fPIDTPCKaon", &fPIDTPCKaon[0], "fPIDTPCKaon[4]/D");
215 fJPsiTree ->Branch("fPIDTPCProton", &fPIDTPCProton[0], "fPIDTPCProton[4]/D");
216
217 fJPsiTree ->Branch("fPIDTOFMuon", &fPIDTOFMuon[0], "fPIDTOFMuon[4]/D");
218 fJPsiTree ->Branch("fPIDTOFElectron", &fPIDTOFElectron[0], "fPIDTOFElectron[4]/D");
219 fJPsiTree ->Branch("fPIDTOFPion", &fPIDTOFPion[0], "fPIDTOFPion[4]/D");
220 fJPsiTree ->Branch("fPIDTOFKaon", &fPIDTOFKaon[0], "fPIDTOFKaon[4]/D");
221 fJPsiTree ->Branch("fPIDTOFProton", &fPIDTOFProton[0], "fPIDTOFProton[4]/D");
8a2486ba 222
489951cd 223 fJPsiTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
224 fJPsiTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
d457ec9e
M		 225 fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
226 fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
227
489951cd 228 fJPsiTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
d0071d7b 229 fJPsiTree ->Branch("fSpdVtxPos", &fSpdVtxPos[0], "fSpdVtxPos[3]/D");
489951cd 230
3d16cd00 231 fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
232 fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
233 fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
234 fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
235 fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
236 fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
237 fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
238 if( fType == 0 ) {
239 fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
240 }
241 if( fType == 1 ) {
242 fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
243 }
8a2486ba 244 if(isMC) {
245 fJPsiTree ->Branch("fGenPart", &fGenPart);
246 fJPsiTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
247 }
248
3d16cd00 249
250 //output tree with Psi2s candidate events
251 fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
252 fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
253 fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
254 fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
255
256 fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
257 fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
258 fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
259 fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
260 fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
489951cd 261 fPsi2sTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
3d16cd00 262 fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
d457ec9e 263
03cff2a0
MB		 264 fPsi2sTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
265 fPsi2sTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
8b17ae4a 266 fPsi2sTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
03cff2a0 267
d9117245 268 fPsi2sTree ->Branch("fPIDTPCMuon", &fPIDTPCMuon[0], "fPIDTPCMuon[4]/D");
269 fPsi2sTree ->Branch("fPIDTPCElectron", &fPIDTPCElectron[0], "fPIDTPCElectron[4]/D");
270 fPsi2sTree ->Branch("fPIDTPCPion", &fPIDTPCPion[0], "fPIDTPCPion[4]/D");
271 fPsi2sTree ->Branch("fPIDTPCKaon", &fPIDTPCKaon[0], "fPIDTPCKaon[4]/D");
272 fPsi2sTree ->Branch("fPIDTPCProton", &fPIDTPCProton[0], "fPIDTPCProton[4]/D");
273
274 fPsi2sTree ->Branch("fPIDTOFMuon", &fPIDTOFMuon[0], "fPIDTOFMuon[4]/D");
275 fPsi2sTree ->Branch("fPIDTOFElectron", &fPIDTOFElectron[0], "fPIDTOFElectron[4]/D");
276 fPsi2sTree ->Branch("fPIDTOFPion", &fPIDTOFPion[0], "fPIDTOFPion[4]/D");
277 fPsi2sTree ->Branch("fPIDTOFKaon", &fPIDTOFKaon[0], "fPIDTOFKaon[4]/D");
278 fPsi2sTree ->Branch("fPIDTOFProton", &fPIDTOFProton[0], "fPIDTOFProton[4]/D");
8a2486ba 279
489951cd 280 fPsi2sTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
281 fPsi2sTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
d457ec9e
M		 282 fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
283 fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
284
489951cd 285 fPsi2sTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
d0071d7b 286 fPsi2sTree ->Branch("fSpdVtxPos", &fSpdVtxPos[0], "fSpdVtxPos[3]/D");
489951cd 287
3d16cd00 288 fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
289 fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
290 fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
291 fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
292 fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
293 fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
294 fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
295 if( fType == 0 ) {
296 fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
297 }
298 if( fType == 1 ) {
299 fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
300 }
8a2486ba 301 if(isMC) {
302 fPsi2sTree ->Branch("fGenPart", &fGenPart);
303 fPsi2sTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
304 }
305
3d16cd00 306
f052ef6f 307 fListTrig = new TList();
308 fListTrig ->SetOwner();
309
4e649ec8 310 fHistCcup4TriggersPerRun = new TH1D("fHistCcup4TriggersPerRun", "fHistCcup4TriggersPerRun", 33000, 167000.5, 200000.5);
311 fListTrig->Add(fHistCcup4TriggersPerRun);
f052ef6f 312
4e649ec8 313 fHistCcup7TriggersPerRun = new TH1D("fHistCcup7TriggersPerRun", "fHistCcup7TriggersPerRun", 33000, 167000.5, 200000.5);
314 fListTrig->Add(fHistCcup7TriggersPerRun);
315
316 fHistCcup2TriggersPerRun = new TH1D("fHistCcup2TriggersPerRun", "fHistCcup2TriggersPerRun", 33000, 167000.5, 200000.5);
317 fListTrig->Add(fHistCcup2TriggersPerRun);
318
adf61188 319 fHistCint1TriggersPerRun = new TH1D("fHistCint1TriggersPerRun", "fHistCint1TriggersPerRun", 33000, 167000.5, 200000.5);
320 fListTrig->Add(fHistCint1TriggersPerRun);
321
a4c57c19 322 fHistC0tvxAndCint1TriggersPerRun = new TH1D("fHistC0tvxAndCint1TriggersPerRun", "fHistC0tvxAndCint1TriggersPerRun", 33000, 167000.5, 200000.5);
323 fListTrig->Add(fHistC0tvxAndCint1TriggersPerRun);
324
4e649ec8 325 fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 33000, 167000.5, 200000.5);
f052ef6f 326 fListTrig->Add(fHistZedTriggersPerRun);
327
4e649ec8 328 fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 33000, 167000.5, 200000.5);
f052ef6f 329 fListTrig->Add(fHistCvlnTriggersPerRun);
330
4e649ec8 331 fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 33000, 167000.5, 200000.5);
e63dc644
MB		 332 fListTrig->Add(fHistMBTriggersPerRun);
333
4e649ec8 334 fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 33000, 167000.5, 200000.5);
e63dc644
MB		 335 fListTrig->Add(fHistCentralTriggersPerRun);
336
4e649ec8 337 fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 33000, 167000.5, 200000.5);
e63dc644
MB		 338 fListTrig->Add(fHistSemiCentralTriggersPerRun);
339
46e1d1dc 340 fListHist = new TList();
341 fListHist ->SetOwner();
342
d1fa9007 343 TString CutNameJPsi[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Two good tracks",
17c65770 344 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
d1fa9007 345 fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",13,0.5,13.5);
346 for (Int_t i = 0; i<13; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
46e1d1dc 347 fListHist->Add(fHistNeventsJPsi);
17c65770 348
349 fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
46e1d1dc 350 fListHist->Add(fHistTPCsignalJPsi);
17c65770 351
352 fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
46e1d1dc 353 fListHist->Add(fHistDiLeptonPtJPsi);
17c65770 354
355 fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
356 fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
46e1d1dc 357 fListHist->Add(fHistDiElectronMass);
17c65770 358
359 fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
360 fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
46e1d1dc 361 fListHist->Add(fHistDiMuonMass);
80b87448 362
363 fHistDiLeptonMass = new TH1D("fHistDiLeptonMass","Invariant mass of J/#psi candidates",130,2.1,6.0);
364 fHistDiLeptonMass->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}) (GeV/c)");
365 fListHist->Add(fHistDiLeptonMass);
17c65770 366
d1fa9007 367 TString CutNamePsi2s[14] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Four good tracks",
17c65770 368 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
369
d1fa9007 370 fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",14,0.5,14.5);
371 for (Int_t i = 0; i<14; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
46e1d1dc 372 fListHist->Add(fHistNeventsPsi2s);
17c65770 373
374 fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
375 fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
376 fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
46e1d1dc 377 fListHist->Add(fHistPsi2sMassVsPt);
17c65770 378
22b569e5 379 fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",50,2.5,5.5);
17c65770 380 fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
46e1d1dc 381 fListHist->Add(fHistPsi2sMassCoherent);
17c65770 382
80b87448 383 fListSystematics = new TList();
384 fListSystematics->SetOwner();
385 fListSystematics->SetName("fListSystematics");
386 fListHist->Add(fListSystematics);
387 InitSystematics();
388
389
3d16cd00 390 PostData(1, fJPsiTree);
391 PostData(2, fPsi2sTree);
f052ef6f 392 PostData(3, fListTrig);
46e1d1dc 393 PostData(4, fListHist);
3d16cd00 394
395}//UserCreateOutputObjects
396
80b87448 397//_____________________________________________________________________________
398void AliAnalysisTaskUpcPsi2s::InitSystematics()
399{
400
401fListJPsiLoose = new TList();
402fListJPsiLoose->SetOwner();
403fListJPsiLoose->SetName("JPsiLoose");
404fListSystematics->Add(fListJPsiLoose);
405
406TH1D *fHistJPsiNClusLoose = new TH1D("JPsiNClusLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
407fListJPsiLoose->Add(fHistJPsiNClusLoose);
408
409TH1D *fHistJPsiChi2Loose = new TH1D("JPsiChi2Loose","Invariant mass of J/#psi candidates",130,2.1,6.0);
410fListJPsiLoose->Add(fHistJPsiChi2Loose);
411
412TH1D *fHistJPsiDCAzLoose = new TH1D("JPsiDCAzLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
413fListJPsiLoose->Add(fHistJPsiDCAzLoose);
414
415TH1D *fHistJPsiDCAxyLoose = new TH1D("JPsiDCAxyLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
416fListJPsiLoose->Add(fHistJPsiDCAxyLoose);
417
22b569e5 418TH1D *fHistJPsiITShitsLoose = new TH1D("JPsiITShitsLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
419fListJPsiLoose->Add(fHistJPsiITShitsLoose);
420
421
80b87448 422fListJPsiTight = new TList();
423fListJPsiTight->SetOwner();
424fListJPsiTight->SetName("JPsiTight");
425fListSystematics->Add(fListJPsiTight);
426
427TH1D *fHistJPsiNClusTight = new TH1D("JPsiNClusTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
428fListJPsiTight->Add(fHistJPsiNClusTight);
429
430TH1D *fHistJPsiChi2Tight = new TH1D("JPsiChi2Tight","Invariant mass of J/#psi candidates",130,2.1,6.0);
431fListJPsiTight->Add(fHistJPsiChi2Tight);
432
433TH1D *fHistJPsiDCAzTight = new TH1D("JPsiDCAzTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
434fListJPsiTight->Add(fHistJPsiDCAzTight);
435
436TH1D *fHistJPsiDCAxyTight = new TH1D("JPsiDCAxyTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
437fListJPsiTight->Add(fHistJPsiDCAxyTight);
438
439
22b569e5 440fListPsi2sLoose = new TList();
441fListPsi2sLoose->SetOwner();
442fListPsi2sLoose->SetName("Psi2sLoose");
443fListSystematics->Add(fListPsi2sLoose);
444
445TH1D *fHistPsi2sNClusLoose = new TH1D("Psi2sNClusLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
446fListPsi2sLoose->Add(fHistPsi2sNClusLoose);
447
448TH1D *fHistPsi2sChi2Loose = new TH1D("Psi2sChi2Loose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
449fListPsi2sLoose->Add(fHistPsi2sChi2Loose);
450
451TH1D *fHistPsi2sDCAzLoose = new TH1D("Psi2sDCAzLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
452fListPsi2sLoose->Add(fHistPsi2sDCAzLoose);
453
454TH1D *fHistPsi2sDCAxyLoose = new TH1D("Psi2sDCAxyLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
455fListPsi2sLoose->Add(fHistPsi2sDCAxyLoose);
456
457TH1D *fHistPsi2sITShitsLoose = new TH1D("Psi2sITShitsLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
458fListPsi2sLoose->Add(fHistPsi2sITShitsLoose);
459
460
461fListPsi2sTight = new TList();
462fListPsi2sTight->SetOwner();
463fListPsi2sTight->SetName("Psi2sTight");
464fListSystematics->Add(fListPsi2sTight);
465
466TH1D *fHistPsi2sNClusTight = new TH1D("Psi2sNClusTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
467fListPsi2sTight->Add(fHistPsi2sNClusTight);
468
469TH1D *fHistPsi2sChi2Tight = new TH1D("Psi2sChi2Tight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
470fListPsi2sTight->Add(fHistPsi2sChi2Tight);
471
472TH1D *fHistPsi2sDCAzTight = new TH1D("Psi2sDCAzTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
473fListPsi2sTight->Add(fHistPsi2sDCAzTight);
474
475TH1D *fHistPsi2sDCAxyTight = new TH1D("Psi2sDCAxyTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
476fListPsi2sTight->Add(fHistPsi2sDCAxyTight);
477
478
80b87448 479}
480
3d16cd00 481//_____________________________________________________________________________
482void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
483{
484
485 //cout<<"#################### Next event ##################"<<endl;
486
f052ef6f 487 if( fType == 0 ){
488 RunESDtrig();
489 if(fRunHist) RunESDhist();
490 if(fRunTree) RunESDtree();
491 }
492
493 if( fType == 1 ){
494 RunAODtrig();
17c65770 495 if(fRunHist) RunAODhist();
496 if(fRunTree) RunAODtree();
497 }
3d16cd00 498
499}//UserExec
17c65770 500//_____________________________________________________________________________
f052ef6f 501void AliAnalysisTaskUpcPsi2s::RunAODtrig()
502{
503
504 //input event
505 AliAODEvent *aod = (AliAODEvent*) InputEvent();
506 if(!aod) return;
507
508 fRunNum = aod ->GetRunNumber();
509 //Trigger
510 TString trigger = aod->GetFiredTriggerClasses();
511
4e649ec8 512 if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
513 if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
514 if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
f052ef6f 515
adf61188 516 if(trigger.Contains("CINT1")) fHistCint1TriggersPerRun->Fill(fRunNum); //CINT1 triggers
a793bb82 517
518 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
519 if(trigger.Contains("CINT1") && (fL0inputs & (1 << 3))) fHistC0tvxAndCint1TriggersPerRun->Fill(fRunNum); //0TVX triggers in CINT1 events
adf61188 520
e63dc644
MB		 521 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
522 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
f052ef6f 523
524 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
525 if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
e63dc644
MB		 526
527 //MB, Central and SemiCentral triggers
528 AliCentrality *centrality = aod->GetCentrality();
529 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
530
03cff2a0
MB		 531 Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
532 //Double_t percentile = centrality->GetCentralityPercentile("V0M");
e63dc644 533
d1fa9007 534 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
535
4e649ec8 536 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
d1fa9007 537
538 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
3c79e118 539
f052ef6f 540PostData(3, fListTrig);
541
542}
543//_____________________________________________________________________________
17c65770 544void AliAnalysisTaskUpcPsi2s::RunAODhist()
545{
546
547 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
548
549 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
550 Double_t muonMass = partMuon->Mass();
551
552 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
553 Double_t electronMass = partElectron->Mass();
554
555 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
556 Double_t pionMass = partPion->Mass();
557
558 //input event
559 AliAODEvent *aod = (AliAODEvent*) InputEvent();
560 if(!aod) return;
80b87448 561
5d08abf3 562 // cout<<"Event number: "<<((TTree*) GetInputData(0))->GetTree()->GetReadEntry()<<endl;
17c65770 563
564 fHistNeventsJPsi->Fill(1);
565 fHistNeventsPsi2s->Fill(1);
566
567 //Trigger
568 TString trigger = aod->GetFiredTriggerClasses();
569
8a2486ba 570 if(!isMC && !trigger.Contains("CCUP") ) return;
17c65770 571
572 fHistNeventsJPsi->Fill(2);
573 fHistNeventsPsi2s->Fill(2);
574
575 //primary vertex
576 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
577 fVtxContrib = fAODVertex->GetNContributors();
578 if(fVtxContrib < 2) return;
579
580 fHistNeventsJPsi->Fill(3);
581 fHistNeventsPsi2s->Fill(3);
582
17c65770 583 //VZERO, ZDC
584 AliAODVZERO *fV0data = aod ->GetVZEROData();
f052ef6f 585 AliAODZDC *fZDCdata = aod->GetZDCData();
17c65770 586
587 fV0Adecision = fV0data->GetV0ADecision();
588 fV0Cdecision = fV0data->GetV0CDecision();
589 if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
590
d1fa9007 591 fHistNeventsJPsi->Fill(4);
592 fHistNeventsPsi2s->Fill(4);
593
f052ef6f 594 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
595 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
596
e63dc644 597 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
f052ef6f 598
d1fa9007 599 fHistNeventsJPsi->Fill(5);
80b87448 600 fHistNeventsPsi2s->Fill(5);
601
602 //Systematics - cut variation
603 if(fRunSystematics) RunAODsystematics(aod);
17c65770 604
17c65770 605 //Two tracks loop
fc8ccb59 606 Int_t nGoodTracks = 0;
17c65770 607 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
608
609 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
610 Short_t qLepton[4], qPion[4];
22b569e5 611 UInt_t nLepton=0, nPion=0, nHighPt=0, nSpdHits=0;
80b87448 612 Double_t fRecTPCsignal[5], fRecTPCsignalDist;
613 Int_t fChannel = 0;
17c65770 614 Int_t mass[3]={-1,-1,-1};
22b569e5 615 Double_t TrackPt[5]={0,0,0,0,0};
616 Double_t MeanPt = -1;
73c6b208 617
fc8ccb59 618
73c6b208 619 //Four track loop
620 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 621 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
73c6b208 622 if( !trk ) continue;
d1fa9007 623 if(!(trk->TestFilterBit(1<<0))) continue;
73c6b208 624
625 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
626 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
627 if(trk->GetTPCNcls() < 50)continue;
628 if(trk->Chi2perNDF() > 4)continue;
629 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
fc8ccb59
MB		 630 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
631 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
632 delete trk_clone;
73c6b208 633 if(TMath::Abs(dca[1]) > 2) continue;
22b569e5 634 Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
635 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
636 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
73c6b208 637
638 TrackIndex[nGoodTracks] = itr;
22b569e5 639 TrackPt[nGoodTracks] = trk->Pt();
73c6b208 640 nGoodTracks++;
641
642 if(nGoodTracks > 4) break;
643 }//Track loop
644
17c65770 645 nLepton=0; nPion=0; nHighPt=0;
646 mass[0]= -1; mass[1]= -1, mass[2]= -1;
647
22b569e5 648 if(nGoodTracks == 4 && nSpdHits>1){
649 MeanPt = GetMedian(TrackPt);
d1fa9007 650 fHistNeventsPsi2s->Fill(6);
17c65770 651 for(Int_t i=0; i<4; i++){
f15c1f69 652 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[i]));
653 if(!trk) AliFatal("Not a standard AOD");
654
22b569e5 655 if(trk->Pt() > MeanPt){
489951cd 656 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
17c65770 657 qLepton[nLepton] = trk->Charge();
489951cd 658 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
17c65770 659 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
660 mass[nLepton] = 0;
661 }
489951cd 662 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
17c65770 663 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
664 mass[nLepton] = 1;
665 }
666 nLepton++;
667 }
668 else{
669 qPion[nPion] = trk->Charge();
670 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
671 nPion++;
672 }
673
674 if(nLepton > 2 || nPion > 2) break;
675 }
676 if((nLepton == 2) && (nPion == 2)){
d1fa9007 677 fHistNeventsPsi2s->Fill(7);
678 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(8);
679 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(9);
680 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
17c65770 681 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
d1fa9007 682 fHistNeventsPsi2s->Fill(11);
22b569e5 683 if(mass[0] != -1 && mass[1] != -1) {
d1fa9007 684 fHistNeventsPsi2s->Fill(12);
17c65770 685 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
686 vDilepton = vLepton[0]+vLepton[1];
687 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
22b569e5 688 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
689 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
690 else {
691 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
692 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
693 }
694
695 if(fChannel == -1) {
696 fHistNeventsPsi2s->Fill(13);
697 if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
698 }
699 if(fChannel == 1){
700 fHistNeventsPsi2s->Fill(14);
701 if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) fHistPsi2sMassCoherent->Fill(vCandidate.M());
702 }
17c65770 703 }
704 }
705 }
706 }
707
fc8ccb59
MB		 708 nGoodTracks = 0;
709 //Two track loop
710 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 711 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
fc8ccb59 712 if( !trk ) continue;
d1fa9007 713 if(!(trk->TestFilterBit(1<<0))) continue;
fc8ccb59
MB		 714
715 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
716 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
717 if(trk->GetTPCNcls() < 70)continue;
718 if(trk->Chi2perNDF() > 4)continue;
719 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
720 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
721 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
722 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
723 delete trk_clone;
724 if(TMath::Abs(dca[1]) > 2) continue;
80b87448 725 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
726 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
fc8ccb59
MB		 727
728 TrackIndex[nGoodTracks] = itr;
729 nGoodTracks++;
730
731 if(nGoodTracks > 2) break;
732 }//Track loop
733
734 nLepton=0; nPion=0; nHighPt=0;
735 mass[0]= -1; mass[1]= -1, mass[2]= -1;
736
737 if(nGoodTracks == 2){
d1fa9007 738 fHistNeventsJPsi->Fill(6);
fc8ccb59 739 for(Int_t i=0; i<2; i++){
f15c1f69 740 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[i]));
741 if(!trk) AliFatal("Not a standard AOD");
fc8ccb59 742 if(trk->Pt() > 1) nHighPt++;
489951cd 743 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
fc8ccb59 744 qLepton[nLepton] = trk->Charge();
489951cd 745 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
fc8ccb59
MB		 746 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
747 mass[nLepton] = 0;
748 }
489951cd 749 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
fc8ccb59
MB		 750 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
751 mass[nLepton] = 1;
752 }
753 nLepton++;
754 }
755 if(nLepton == 2){
d1fa9007 756 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(7);
fc8ccb59 757 if(qLepton[0]*qLepton[1] < 0){
d1fa9007 758 fHistNeventsJPsi->Fill(8);
fc8ccb59 759 if(nHighPt > 0){
d1fa9007 760 fHistNeventsJPsi->Fill(9);
489951cd 761 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
d1fa9007 762 if(nHighPt == 2) fHistNeventsJPsi->Fill(10);
80b87448 763 if(mass[0] != -1 && mass[1] != -1) {
d1fa9007 764 fHistNeventsJPsi->Fill(11);
fc8ccb59 765 vCandidate = vLepton[0]+vLepton[1];
80b87448 766 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
767 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
768 else {
769 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
770 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
771 }
fc8ccb59 772 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
80b87448 773 if(fChannel == -1) {
fc8ccb59 774 fHistDiMuonMass->Fill(vCandidate.M());
80b87448 775 if(vCandidate.Pt()<0.15)fHistDiLeptonMass->Fill(vCandidate.M());
d1fa9007 776 fHistNeventsJPsi->Fill(12);
fc8ccb59 777 }
80b87448 778 if(fChannel == 1) {
fc8ccb59 779 fHistDiElectronMass->Fill(vCandidate.M());
80b87448 780 if(vCandidate.Pt()<0.3)fHistDiLeptonMass->Fill(vCandidate.M());
d1fa9007 781 fHistNeventsJPsi->Fill(13);
fc8ccb59
MB		 782 }
783 }
784 }
785 }
786 }
787 }
788
789
46e1d1dc 790 PostData(4, fListHist);
17c65770 791
792}
3d16cd00 793
80b87448 794//_____________________________________________________________________________
d9117245 795void AliAnalysisTaskUpcPsi2s::RunAODtree()
80b87448 796{
d9117245 797 //input event
798 AliAODEvent *aod = (AliAODEvent*) InputEvent();
799 if(!aod) return;
80b87448 800
d9117245 801 if(isMC) RunAODMC(aod);
80b87448 802
d9117245 803 //input data
804 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
805 fDataFilnam->Clear();
806 fDataFilnam->SetString(filnam);
807 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
808 fRunNum = aod ->GetRunNumber();
80b87448 809
d9117245 810 //Trigger
811 TString trigger = aod->GetFiredTriggerClasses();
80b87448 812
d9117245 813 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
814 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
815 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
adf61188 816 fTrigger[3] = trigger.Contains("CINT1");
d9117245 817
818 Bool_t isTriggered = kFALSE;
819 for(Int_t i=0; i<ntrg; i++) {
820 if( fTrigger[i] ) isTriggered = kTRUE;
821 }
822 if(!isMC && !isTriggered ) return;
80b87448 823
d9117245 824 //trigger inputs
825 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
826 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
80b87448 827
d9117245 828 //Event identification
829 fPerNum = aod ->GetPeriodNumber();
830 fOrbNum = aod ->GetOrbitNumber();
831 fBCrossNum = aod ->GetBunchCrossNumber();
80b87448 832
d9117245 833 //primary vertex
834 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
835 fVtxContrib = fAODVertex->GetNContributors();
836 fVtxPos[0] = fAODVertex->GetX();
837 fVtxPos[1] = fAODVertex->GetY();
838 fVtxPos[2] = fAODVertex->GetZ();
839 Double_t CovMatx[6];
840 fAODVertex->GetCovarianceMatrix(CovMatx);
841 fVtxErr[0] = CovMatx[0];
842 fVtxErr[1] = CovMatx[1];
843 fVtxErr[2] = CovMatx[2];
844 fVtxChi2 = fAODVertex->GetChi2();
845 fVtxNDF = fAODVertex->GetNDF();
d0071d7b 846
847 //SPD primary vertex
848 AliAODVertex *fSPDVertex = aod->GetPrimaryVertexSPD();
dd9393a2 849 if(fSPDVertex->GetNContributors()>0){
d0071d7b 850 fSpdVtxPos[0] = fSPDVertex->GetX();
851 fSpdVtxPos[1] = fSPDVertex->GetY();
852 fSpdVtxPos[2] = fSPDVertex->GetZ();
853 }
d4d2400a 854 else{
855 fSpdVtxPos[0] = -666;
856 fSpdVtxPos[1] = -666;
857 fSpdVtxPos[2] = -666;
858 }
80b87448 859
d9117245 860 //Tracklets
861 fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
80b87448 862
d9117245 863 //VZERO, ZDC
864 AliAODVZERO *fV0data = aod ->GetVZEROData();
865 AliAODZDC *fZDCdata = aod->GetZDCData();
80b87448 866
d9117245 867 fV0Adecision = fV0data->GetV0ADecision();
868 fV0Cdecision = fV0data->GetV0CDecision();
869 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
870 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
80b87448 871
d9117245 872 fNLooseTracks = 0;
80b87448 873
d9117245 874 //Track loop - loose cuts
875 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 876 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
d9117245 877 if( !trk ) continue;
878 if(!(trk->TestFilterBit(1<<0))) continue;
80b87448 879
d9117245 880 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
881 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
882 if(trk->GetTPCNcls() < 20)continue;
883 fNLooseTracks++;
884 }//Track loop -loose cuts
885
886 Int_t nGoodTracks=0;
887 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
80b87448 888
80b87448 889 //Two track loop
890 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 891 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
80b87448 892 if( !trk ) continue;
893 if(!(trk->TestFilterBit(1<<0))) continue;
d9117245 894
895 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
896 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
897 if(trk->GetTPCNcls() < 70)continue;
898 if(trk->Chi2perNDF() > 4)continue;
899 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
80b87448 900 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
901 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
902 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
903 delete trk_clone;
d9117245 904 if(TMath::Abs(dca[1]) > 2) continue;
80b87448 905 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
d9117245 906 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
80b87448 907
80b87448 908
909 TrackIndex[nGoodTracks] = itr;
910 nGoodTracks++;
911
912 if(nGoodTracks > 2) break;
913 }//Track loop
80b87448 914
d9117245 915 fJPsiAODTracks->Clear("C");
80b87448 916 if(nGoodTracks == 2){
d9117245 917
918 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
919 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
920 Double_t muonMass = partMuon->Mass();
921 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
922 Double_t electronMass = partElectron->Mass();
923 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
924 Double_t pionMass = partPion->Mass();
925
926 Double_t KFcov[21];
927 Double_t KFpar[6];
928 Double_t KFmass = pionMass;
929 Double_t fRecTPCsignal;
930 AliKFParticle *KFpart[2];
931 AliKFVertex *KFvtx = new AliKFVertex();
932 KFvtx->SetField(aod->GetMagneticField());
933
934 for(Int_t i=0; i<2; i++){
f15c1f69 935 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[i]));
936 if(!trk) AliFatal("Not a standard AOD");
937
d9117245 938 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
939 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
940 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
941 delete trk_clone;
942
943 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
944 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
945
946 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
947 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
948 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
949 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
950 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
951
952 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
953 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
954 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
955 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
956 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
957
958 trk->GetPosition(KFpar);
959 trk->PxPyPz(KFpar+3);
960 trk->GetCovMatrix(KFcov);
961
962 if(trk->Pt() > 1){
963 fRecTPCsignal = trk->GetTPCsignal();
964 if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
965 if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
80b87448 966 }
d9117245 967 else KFmass = pionMass;
968
969 KFpart[i] = new AliKFParticle();
970 KFpart[i]->SetField(aod->GetMagneticField());
971 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
972 KFvtx->AddDaughter(*KFpart[i]);
973
974
975 Double_t pos[3]={0,0,0};
976 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
977 parTrk->CopyFromVTrack((AliVTrack*) trk);
978 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
979 else {
980 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
981 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
982 }
983 delete parTrk;
984 }
985 fKfVtxPos[0]= KFvtx->GetX();
986 fKfVtxPos[1]= KFvtx->GetY();
987 fKfVtxPos[2]= KFvtx->GetZ();
988 for(UInt_t i=0; i<2; i++)delete KFpart[i];
989 delete KFvtx;
80b87448 990
d9117245 991 if(!isMC) fJPsiTree ->Fill();
992 }
993
994 nGoodTracks = 0;
995 //Four track loop
80b87448 996 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 997 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
80b87448 998 if( !trk ) continue;
999 if(!(trk->TestFilterBit(1<<0))) continue;
1000
d9117245 1001 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
1002 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
1003 if(trk->GetTPCNcls() < 50)continue;
1004 if(trk->Chi2perNDF() > 4)continue;
80b87448 1005 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1006 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1007 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1008 delete trk_clone;
d9117245 1009 if(TMath::Abs(dca[1]) > 2) continue;
1010 Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1011 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
1012
80b87448 1013 TrackIndex[nGoodTracks] = itr;
1014 nGoodTracks++;
1015
d9117245 1016 if(nGoodTracks > 4) break;
80b87448 1017 }//Track loop
d9117245 1018
1019 fPsi2sAODTracks->Clear("C");
1020 if(nGoodTracks == 4){
1021
1022 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1023 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1024 Double_t muonMass = partMuon->Mass();
1025 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1026 Double_t electronMass = partElectron->Mass();
1027 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1028 Double_t pionMass = partPion->Mass();
1029
1030 Double_t KFcov[21];
1031 Double_t KFpar[6];
1032 Double_t KFmass = pionMass;
1033 Double_t fRecTPCsignal;
1034 AliKFParticle *KFpart[4];
1035 AliKFVertex *KFvtx = new AliKFVertex();
1036 KFvtx->SetField(aod->GetMagneticField());
1037
1038 for(Int_t i=0; i<4; i++){
f15c1f69 1039 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[i]));
1040 if(!trk) AliFatal("Not a standard AOD");
1041
d9117245 1042 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1043 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1044 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1045 delete trk_clone;
1046
1047 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
1048 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
1049
1050
1051 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1052 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1053 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1054 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1055 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
1056
1057 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1058 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1059 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1060 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1061 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
1062
1063 trk->GetPosition(KFpar);
1064 trk->PxPyPz(KFpar+3);
1065 trk->GetCovMatrix(KFcov);
1066
1067 if(trk->Pt() > 1){
1068 fRecTPCsignal = trk->GetTPCsignal();
1069 if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
1070 if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
1071 }
1072 else KFmass = pionMass;
1073
1074 KFpart[i] = new AliKFParticle();
1075 KFpart[i]->SetField(aod->GetMagneticField());
1076 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
1077 KFvtx->AddDaughter(*KFpart[i]);
1078
1079 Double_t pos[3]={0,0,0};
1080 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
1081 parTrk->CopyFromVTrack((AliVTrack*) trk);
1082 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
1083 else {
1084 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1085 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1086 }
1087 delete parTrk;
1088 }
1089 fKfVtxPos[0]= KFvtx->GetX();
1090 fKfVtxPos[1]= KFvtx->GetY();
1091 fKfVtxPos[2]= KFvtx->GetZ();
1092 for(UInt_t i=0; i<4; i++)delete KFpart[i];
1093 delete KFvtx;
1094 if(!isMC) fPsi2sTree ->Fill();
1095 }
1096
1097 if(isMC){
1098 fJPsiTree ->Fill();
1099 fPsi2sTree ->Fill();
1100 }
1101
1102 PostData(1, fJPsiTree);
1103 PostData(2, fPsi2sTree);
1104
1105}//RunAOD
1106
1107
1108//_____________________________________________________________________________
1109void AliAnalysisTaskUpcPsi2s::RunAODMC(AliAODEvent *aod)
1110{
1111
c775e261 1112 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
1113 fTriggerInputsMC[0] = kFALSE;//0SM2
1114 fTriggerInputsMC[1] = fL0inputs & (1 << 0);//0VBA
1115 fTriggerInputsMC[2] = fL0inputs & (1 << 1);//0VBC
1116 fTriggerInputsMC[3] = fL0inputs & (1 << 9);//0OMU
1117
d9117245 1118 fGenPart->Clear("C");
1119
1120 TClonesArray *arrayMC = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1121 if(!arrayMC) return;
1122
1123 Int_t nmc=0;
1124 //loop over mc particles
1125 for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
1126 AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(imc);
1127 if(!mcPart) continue;
1128
1129 if(mcPart->GetMother() >= 0) continue;
1130
1131 TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1132 part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1133 part->SetPdgCode(mcPart->GetPdgCode());
1134 part->SetUniqueID(imc);
1135 }//loop over mc particles
1136
1137}//RunAODMC
1138
1139
1140//_____________________________________________________________________________
1141void AliAnalysisTaskUpcPsi2s::RunESDtrig()
1142{
1143
1144 //input event
1145 AliESDEvent *esd = (AliESDEvent*) InputEvent();
1146 if(!esd) return;
1147
1148 fRunNum = esd ->GetRunNumber();
1149 //Trigger
1150 TString trigger = esd->GetFiredTriggerClasses();
1151
1152 if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
1153 if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
1154 if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
1155
a4c57c19 1156 if(trigger.Contains("CINT1")) fHistCint1TriggersPerRun->Fill(fRunNum); //CINT1 triggers
1157 if(trigger.Contains("CINT1") && trigger.Contains("C0TVX")) fHistC0tvxAndCint1TriggersPerRun->Fill(fRunNum); //C0TVX triggers in CINT1 events
1158
d9117245 1159 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
1160 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
1161
1162 if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
1163
1164 //MB, Central and SemiCentral triggers
1165 AliCentrality *centrality = esd->GetCentrality();
1166 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1167
1168 //Double_t percentile = centrality->GetCentralityPercentile("V0M");
1169 Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
1170
1171 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
1172
1173 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
1174
1175 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
1176
1177
1178PostData(3, fListTrig);
1179
1180}
1181//_____________________________________________________________________________
1182void AliAnalysisTaskUpcPsi2s::RunESDhist()
1183{
1184
1185 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1186
1187 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1188 Double_t muonMass = partMuon->Mass();
1189
1190 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1191 Double_t electronMass = partElectron->Mass();
1192
1193 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1194 Double_t pionMass = partPion->Mass();
1195
1196 //input event
1197 AliESDEvent *esd = (AliESDEvent*) InputEvent();
1198 if(!esd) return;
1199
1200 fHistNeventsJPsi->Fill(1);
1201 fHistNeventsPsi2s->Fill(1);
1202
1203 //Trigger
1204 TString trigger = esd->GetFiredTriggerClasses();
1205
1206 if(!isMC && !trigger.Contains("CCUP") ) return;
1207
1208 fHistNeventsJPsi->Fill(2);
1209 fHistNeventsPsi2s->Fill(2);
1210
1211 //primary vertex
1212 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1213 fVtxContrib = fESDVertex->GetNContributors();
1214 if(fVtxContrib < 2) return;
1215
1216 fHistNeventsJPsi->Fill(3);
1217 fHistNeventsPsi2s->Fill(3);
1218
1219 //VZERO, ZDC
1220 AliESDVZERO *fV0data = esd->GetVZEROData();
1221 AliESDZDC *fZDCdata = esd->GetESDZDC();
1222
1223 fV0Adecision = fV0data->GetV0ADecision();
1224 fV0Cdecision = fV0data->GetV0CDecision();
1225 if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
1226
1227 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1228 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1229 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
1230
1231 fHistNeventsJPsi->Fill(4);
1232 fHistNeventsPsi2s->Fill(4);
1233
1234 Int_t nGoodTracks=0;
1235 //Two tracks loop
1236 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1237
1238 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1239 Short_t qLepton[4], qPion[4];
1240 UInt_t nLepton=0, nPion=0, nHighPt=0;
1241 Double_t fRecTPCsignal[5];
80b87448 1242 Int_t mass[3]={-1,-1,-1};
d9117245 1243
1244 //Two Track loop
1245 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1246 AliESDtrack *trk = esd->GetTrack(itr);
1247 if( !trk ) continue;
1248
1249 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1250 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1251 if(trk->GetTPCNcls() < 70)continue;
1252 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1253 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1254 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1255 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1256 trk->GetImpactParameters(dca[0],dca[1]);
1257 if(TMath::Abs(dca[1]) > 2) continue;
1258 if(TMath::Abs(dca[1]) > 0.2) continue;
1259
1260 TrackIndex[nGoodTracks] = itr;
1261 nGoodTracks++;
1262 if(nGoodTracks > 2) break;
1263 }//Track loop
1264
80b87448 1265
1266 if(nGoodTracks == 2){
d9117245 1267 fHistNeventsJPsi->Fill(5);
1268 for(Int_t i=0; i<2; i++){
1269 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
80b87448 1270 if(trk->Pt() > 1) nHighPt++;
1271 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1272 qLepton[nLepton] = trk->Charge();
1273 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1274 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1275 mass[nLepton] = 0;
1276 }
1277 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1278 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1279 mass[nLepton] = 1;
1280 }
1281 nLepton++;
1282 }
1283 if(nLepton == 2){
d9117245 1284 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
1285 if(qLepton[0]*qLepton[1] < 0){
1286 fHistNeventsJPsi->Fill(7);
1287 if(nHighPt > 0){
1288 fHistNeventsJPsi->Fill(8);
1289 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
1290 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
1291 if(mass[0] == mass[1] && mass[0] != -1) {
1292 fHistNeventsJPsi->Fill(10);
1293 vCandidate = vLepton[0]+vLepton[1];
1294 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
1295 if(mass[0] == 0) {
1296 fHistDiMuonMass->Fill(vCandidate.M());
1297 fHistNeventsJPsi->Fill(11);
1298 }
1299 if(mass[0] == 1) {
1300 fHistDiElectronMass->Fill(vCandidate.M());
1301 fHistNeventsJPsi->Fill(12);
1302 }
22b569e5 1303 }
d9117245 1304 }
22b569e5 1305 }
d9117245 1306 }
1307 }
1308 nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1309 mass[0]= -1; mass[1]= -1, mass[2]= -1;
1310
1311 //Four Track loop
1312 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1313 AliESDtrack *trk = esd->GetTrack(itr);
22b569e5 1314 if( !trk ) continue;
22b569e5 1315
1316 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1317 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
d9117245 1318 if(trk->GetTPCNcls() < 50)continue;
1319 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1320 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1321 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1322 trk->GetImpactParameters(dca[0],dca[1]);
1323 if(TMath::Abs(dca[1]) > 2) continue;
22b569e5 1324
22b569e5 1325 TrackIndex[nGoodTracks] = itr;
22b569e5 1326 nGoodTracks++;
d9117245 1327 if(nGoodTracks > 4) break;
22b569e5 1328 }//Track loop
22b569e5 1329
1330 if(nGoodTracks == 4){
d9117245 1331 fHistNeventsPsi2s->Fill(5);
1332 for(Int_t i=0; i<4; i++){
1333 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1334
1335 if(trk->Pt() > 1){
22b569e5 1336 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1337 qLepton[nLepton] = trk->Charge();
1338 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1339 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1340 mass[nLepton] = 0;
1341 }
1342 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1343 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1344 mass[nLepton] = 1;
1345 }
1346 nLepton++;
1347 }
1348 else{
1349 qPion[nPion] = trk->Charge();
1350 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1351 nPion++;
d9117245 1352 }
1353
1354 if(nLepton > 2 || nPion > 2) break;
22b569e5 1355 }
d9117245 1356 if((nLepton == 2) && (nPion == 2)){
1357 fHistNeventsPsi2s->Fill(6);
1358 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1359 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1360 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1361 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1362 fHistNeventsPsi2s->Fill(10);
1363 if(mass[0] == mass[1]) {
1364 fHistNeventsPsi2s->Fill(11);
1365 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1366 vDilepton = vLepton[0]+vLepton[1];
1367 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1368 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1369 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1370 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1371 }
1372 }
1373 }
1374 }
1375
1376 PostData(4, fListHist);
80b87448 1377
1378}
d9117245 1379
3d16cd00 1380//_____________________________________________________________________________
d9117245 1381void AliAnalysisTaskUpcPsi2s::RunESDtree()
3d16cd00 1382{
3d16cd00 1383
d9117245 1384 //input event
1385 AliESDEvent *esd = (AliESDEvent*) InputEvent();
1386 if(!esd) return;
1387
1388 if(isMC) RunESDMC(esd);
8a2486ba 1389
3d16cd00 1390 //input data
1391 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1392 fDataFilnam->Clear();
1393 fDataFilnam->SetString(filnam);
1394 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
d9117245 1395 fRunNum = esd->GetRunNumber();
3d16cd00 1396
d9117245 1397 //Trigger
1398 TString trigger = esd->GetFiredTriggerClasses();
3d16cd00 1399
8b17ae4a
MB		 1400 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1401 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
1402 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1403
1404 Bool_t isTriggered = kFALSE;
1405 for(Int_t i=0; i<ntrg; i++) {
1406 if( fTrigger[i] ) isTriggered = kTRUE;
1407 }
8a2486ba 1408 if(!isMC && !isTriggered ) return;
d9117245 1409
3d16cd00 1410 //trigger inputs
d9117245 1411 fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1412 fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1413
1414 //TOF trigger info (0OMU)
1415 fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1416 fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
3d16cd00 1417
1418 //Event identification
d9117245 1419 fPerNum = esd->GetPeriodNumber();
1420 fOrbNum = esd->GetOrbitNumber();
1421 fBCrossNum = esd->GetBunchCrossNumber();
3d16cd00 1422
1423 //primary vertex
d9117245 1424 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1425 fVtxContrib = fESDVertex->GetNContributors();
1426 fVtxPos[0] = fESDVertex->GetX();
1427 fVtxPos[1] = fESDVertex->GetY();
1428 fVtxPos[2] = fESDVertex->GetZ();
d457ec9e 1429 Double_t CovMatx[6];
d9117245 1430 fESDVertex->GetCovarianceMatrix(CovMatx);
489951cd 1431 fVtxErr[0] = CovMatx[0];
1432 fVtxErr[1] = CovMatx[1];
1433 fVtxErr[2] = CovMatx[2];
d9117245 1434 fVtxChi2 = fESDVertex->GetChi2();
1435 fVtxNDF = fESDVertex->GetNDF();
d0071d7b 1436
1437 //SPD primary vertex
1438 AliESDVertex *fSPDVertex = (AliESDVertex*) esd->GetPrimaryVertexSPD();
1439 if(fSPDVertex){
1440 fSpdVtxPos[0] = fSPDVertex->GetX();
1441 fSpdVtxPos[1] = fSPDVertex->GetY();
1442 fSpdVtxPos[2] = fSPDVertex->GetZ();
1443 }
3d16cd00 1444
1445 //Tracklets
d9117245 1446 fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
3d16cd00 1447
1448 //VZERO, ZDC
d9117245 1449 AliESDVZERO *fV0data = esd->GetVZEROData();
1450 AliESDZDC *fZDCdata = esd->GetESDZDC();
3d16cd00 1451
1452 fV0Adecision = fV0data->GetV0ADecision();
1453 fV0Cdecision = fV0data->GetV0CDecision();
d9117245 1454 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1455 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1456
489951cd 1457 fNLooseTracks = 0;
1458
1459 //Track loop - loose cuts
d9117245 1460 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1461 AliESDtrack *trk = esd->GetTrack(itr);
489951cd 1462 if( !trk ) continue;
1463
d9117245 1464 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1465 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
489951cd 1466 if(trk->GetTPCNcls() < 20)continue;
1467 fNLooseTracks++;
1468 }//Track loop -loose cuts
1469
17c65770 1470 Int_t nGoodTracks=0;
3d16cd00 1471 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
fc8ccb59 1472
d9117245 1473 //Two Track loop
1474 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1475 AliESDtrack *trk = esd->GetTrack(itr);
3d16cd00 1476 if( !trk ) continue;
d9117245 1477
1478 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1479 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
8b17ae4a 1480 if(trk->GetTPCNcls() < 70)continue;
d9117245 1481 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
8b17ae4a 1482 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
d9117245 1483 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1484 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1485 trk->GetImpactParameters(dca[0],dca[1]);
3d16cd00 1486 if(TMath::Abs(dca[1]) > 2) continue;
d9117245 1487 if(TMath::Abs(dca[0]) > 0.2) continue;
38677cd7 1488
3d16cd00 1489 TrackIndex[nGoodTracks] = itr;
1490 nGoodTracks++;
d9117245 1491 if(nGoodTracks > 2) break;
3d16cd00 1492 }//Track loop
d9117245 1493
1494 fJPsiESDTracks->Clear("C");
8b17ae4a
MB		 1495 if(nGoodTracks == 2){
1496 for(Int_t i=0; i<2; i++){
d9117245 1497 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
3d16cd00 1498
d9117245 1499 AliExternalTrackParam cParam;
1500 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
8b17ae4a 1501
d9117245 1502 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
489951cd 1503
d9117245 1504 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1505 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1506 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1507 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1508 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
489951cd 1509
d9117245 1510 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1511 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1512 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1513 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1514 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
489951cd 1515
8b17ae4a 1516 Double_t pos[3]={0,0,0};
d9117245 1517 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
8b17ae4a
MB		 1518 else {
1519 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1520 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
d9117245 1521 }
3d16cd00 1522 }
8a2486ba 1523 if(!isMC) fJPsiTree ->Fill();
3d16cd00 1524 }
8b17ae4a 1525
d9117245 1526 nGoodTracks = 0;
1527 //Four track loop
1528 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1529 AliESDtrack *trk = esd->GetTrack(itr);
73c6b208 1530 if( !trk ) continue;
1531
d9117245 1532 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1533 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
8b17ae4a 1534 if(trk->GetTPCNcls() < 50)continue;
d9117245 1535 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1536 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1537 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1538 trk->GetImpactParameters(dca[0],dca[1]);
73c6b208 1539 if(TMath::Abs(dca[1]) > 2) continue;
d9117245 1540
73c6b208 1541 TrackIndex[nGoodTracks] = itr;
1542 nGoodTracks++;
d9117245 1543 if(nGoodTracks > 4) break;
73c6b208 1544 }//Track loop
489951cd 1545
d9117245 1546 fPsi2sESDTracks->Clear("C");
1547 if(nGoodTracks == 4){
8b17ae4a 1548 for(Int_t i=0; i<4; i++){
d9117245 1549 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
8a2486ba 1550
d9117245 1551 AliExternalTrackParam cParam;
1552 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1553
1554 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
489951cd 1555
d9117245 1556 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1557 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1558 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1559 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1560 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
489951cd 1561
d9117245 1562 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1563 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1564 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1565 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1566 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
1567
8b17ae4a 1568 Double_t pos[3]={0,0,0};
d9117245 1569 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
8b17ae4a
MB		 1570 else {
1571 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1572 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
d9117245 1573 }
3d16cd00 1574 }
8a2486ba 1575 if(!isMC) fPsi2sTree ->Fill();
1576 }
1577
1578 if(isMC){
1579 fJPsiTree ->Fill();
1580 fPsi2sTree ->Fill();
3d16cd00 1581 }
d9117245 1582
03cff2a0
MB		 1583 PostData(1, fJPsiTree);
1584 PostData(2, fPsi2sTree);
3d16cd00 1585
d9117245 1586}//RunESD
3d16cd00 1587
8a2486ba 1588
1589//_____________________________________________________________________________
d9117245 1590void AliAnalysisTaskUpcPsi2s::RunESDMC(AliESDEvent* esd)
8a2486ba 1591{
1592
d9117245 1593 AliTriggerAnalysis *fTrigAna = new AliTriggerAnalysis();
1594 fTrigAna->SetAnalyzeMC(isMC);
1595
1596 if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) > 1) fTriggerInputsMC[0] = kTRUE;
1597 if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) < 2) fTriggerInputsMC[0] = kFALSE;
1598
1599 fTriggerInputsMC[1] = esd->GetHeader()->IsTriggerInputFired("0OMU");
1600 fTriggerInputsMC[2] = esd->GetHeader()->IsTriggerInputFired("0VBA");
1601 fTriggerInputsMC[3] = esd->GetHeader()->IsTriggerInputFired("0VBC");
1602
8a2486ba 1603 fGenPart->Clear("C");
1604
d9117245 1605 AliMCEvent *mc = MCEvent();
1606 if(!mc) return;
8a2486ba 1607
d9117245 1608 Int_t nmc = 0;
8a2486ba 1609 //loop over mc particles
d9117245 1610 for(Int_t imc=0; imc<mc->GetNumberOfTracks(); imc++) {
1611 AliMCParticle *mcPart = (AliMCParticle*) mc->GetTrack(imc);
8a2486ba 1612 if(!mcPart) continue;
1613
1614 if(mcPart->GetMother() >= 0) continue;
1615
1616 TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1617 part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
d9117245 1618 part->SetPdgCode(mcPart->PdgCode());
8a2486ba 1619 part->SetUniqueID(imc);
1620 }//loop over mc particles
1621
d9117245 1622}//RunESDMC
1623
8a2486ba 1624
1625
3d16cd00 1626//_____________________________________________________________________________
d9117245 1627void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
f052ef6f 1628{
1629
d9117245 1630 cout<<"Analysis complete."<<endl;
1631}//Terminate
03cff2a0 1632
d9117245 1633//_____________________________________________________________________________
1634Double_t AliAnalysisTaskUpcPsi2s::GetMedian(Double_t *daArray) {
1635 // Allocate an array of the same size and sort it.
1636 Double_t dpSorted[4];
1637 for (Int_t i = 0; i < 4; ++i) {
1638 dpSorted[i] = daArray[i];
1639 }
1640 for (Int_t i = 3; i > 0; --i) {
1641 for (Int_t j = 0; j < i; ++j) {
1642 if (dpSorted[j] > dpSorted[j+1]) {
1643 Double_t dTemp = dpSorted[j];
1644 dpSorted[j] = dpSorted[j+1];
1645 dpSorted[j+1] = dTemp;
1646 }
1647 }
1648 }
f052ef6f 1649
d9117245 1650 // Middle or average of middle values in the sorted array.
1651 Double_t dMedian = 0.0;
1652 dMedian = (dpSorted[2] + dpSorted[1])/2.0;
1653
1654 return dMedian;
f052ef6f 1655}
d9117245 1656
f052ef6f 1657//_____________________________________________________________________________
d9117245 1658void AliAnalysisTaskUpcPsi2s::RunAODsystematics(AliAODEvent* aod)
f052ef6f 1659{
1660
d9117245 1661 Double_t fJPsiSels[4];
f052ef6f 1662
d9117245 1663 fJPsiSels[0] = 70; //min number of TPC clusters
1664 fJPsiSels[1] = 4; //chi2
1665 fJPsiSels[2] = 2; //DCAz
1666 fJPsiSels[3] = 1; // DCAxy 1x
f052ef6f 1667
d9117245 1668 Double_t fJPsiSelsMid[4];
f052ef6f 1669
d9117245 1670 fJPsiSelsMid[0] = 70; //min number of TPC clusters
1671 fJPsiSelsMid[1] = 4; //chi2
1672 fJPsiSelsMid[2] = 2; //DCAz
1673 fJPsiSelsMid[3] = 1; // DCAxy 1x
f052ef6f 1674
d9117245 1675 Double_t fJPsiSelsLoose[4];
f052ef6f 1676
d9117245 1677 fJPsiSelsLoose[0] = 60; //min number of TPC clusters
1678 fJPsiSelsLoose[1] = 5; //chi2
1679 fJPsiSelsLoose[2] = 3; //DCAz
1680 fJPsiSelsLoose[3] = 2; // DCAxy 2x
f052ef6f 1681
d9117245 1682 Double_t fJPsiSelsTight[4];
f052ef6f 1683
d9117245 1684 fJPsiSelsTight[0] = 80; //min number of TPC clusters
1685 fJPsiSelsTight[1] = 3.5; //chi2
1686 fJPsiSelsTight[2] = 1; //DCAz
1687 fJPsiSelsTight[3] = 0.5; // DCAxy 0.5x
1688
1689 Int_t nGoodTracks = 0;
f052ef6f 1690 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1691
1692 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
d9117245 1693 Short_t qLepton[4],qPion[4];
f052ef6f 1694 UInt_t nLepton=0, nPion=0, nHighPt=0;
d9117245 1695 Double_t fRecTPCsignal[5], fRecTPCsignalDist;
1696 Int_t fChannel = 0;
f052ef6f 1697
d9117245 1698 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
1699
1700 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1701
1702 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1703 Double_t muonMass = partMuon->Mass();
1704
1705 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1706 Double_t electronMass = partElectron->Mass();
1707
1708 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1709 Double_t pionMass = partPion->Mass();
1710
1711
1712for(Int_t i=0; i<5; i++){
1713 //cout<<"Loose sytematics, cut"<<i<<endl;
1714 for(Int_t j=0; j<4; j++){
1715 if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
1716 else fJPsiSels[j] = fJPsiSelsMid[j];
1717 }
1718 //Two track loop
1719 nGoodTracks = 0;
1720 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 1721 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
f052ef6f 1722 if( !trk ) continue;
d9117245 1723 if(!(trk->TestFilterBit(1<<0))) continue;
f052ef6f 1724
1725 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1726 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
d9117245 1727 if(i!=4){ if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;}
1728 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1729 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1730 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1731 delete trk_clone;
1732 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
f052ef6f 1733
d9117245 1734 if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1735 if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1736 if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1737 if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1738
f052ef6f 1739 TrackIndex[nGoodTracks] = itr;
1740 nGoodTracks++;
d9117245 1741
1742 if(nGoodTracks > 2) break;
f052ef6f 1743 }//Track loop
d9117245 1744
1745 Int_t mass[3]={-1,-1,-1};
1746 fChannel = 0;
1747 nLepton=0; nHighPt=0;
f052ef6f 1748
1749 if(nGoodTracks == 2){
d9117245 1750 for(Int_t k=0; k<2; k++){
f15c1f69 1751 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[k]));
1752 if(!trk) AliFatal("Not a standard AOD");
1753
f052ef6f 1754 if(trk->Pt() > 1) nHighPt++;
489951cd 1755 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
f052ef6f 1756 qLepton[nLepton] = trk->Charge();
489951cd 1757 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
f052ef6f 1758 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1759 mass[nLepton] = 0;
1760 }
489951cd 1761 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
f052ef6f 1762 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1763 mass[nLepton] = 1;
1764 }
1765 nLepton++;
1766 }
1767 if(nLepton == 2){
d9117245 1768 if(qLepton[0]*qLepton[1] < 0 && nHighPt > 0 && (mass[0]!=-1 || mass[1]!=-1)){
1769 vCandidate = vLepton[0]+vLepton[1];
1770 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1771 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1772 else {
1773 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1774 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
f052ef6f 1775 }
d9117245 1776 if(fChannel == -1 && vCandidate.Pt()<0.15) ((TH1D*)(fListJPsiLoose->At(i)))->Fill(vCandidate.M());
1777 if(fChannel == 1 && vCandidate.Pt()<0.3) ((TH1D*)(fListJPsiLoose->At(i)))->Fill(vCandidate.M());
f052ef6f 1778 }
1779 }
1780 }
d9117245 1781}//loose cuts
1782
1783for(Int_t i=0; i<4; i++){
1784 //cout<<"Tight sytematics, cut"<<i<<endl;
1785 for(Int_t j=0; j<4; j++){
1786 if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
1787 else fJPsiSels[j] = fJPsiSelsMid[j];
1788 }
1789 //Two track loop
1790 nGoodTracks = 0;
1791 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 1792 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
03cff2a0 1793 if( !trk ) continue;
d9117245 1794 if(!(trk->TestFilterBit(1<<0))) continue;
03cff2a0
MB		 1795
1796 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1797 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
d9117245 1798 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1799 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1800 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1801 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1802 delete trk_clone;
1803 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
03cff2a0 1804
d9117245 1805 if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1806 if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1807 if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1808 if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1809
03cff2a0
MB		 1810 TrackIndex[nGoodTracks] = itr;
1811 nGoodTracks++;
d9117245 1812
1813 if(nGoodTracks > 2) break;
03cff2a0 1814 }//Track loop
d9117245 1815
1816 Int_t mass[3]={-1,-1,-1};
1817 fChannel = 0;
1818 nLepton=0; nHighPt=0;
1819
1820 if(nGoodTracks == 2){
1821 for(Int_t k=0; k<2; k++){
f15c1f69 1822 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[k]));
1823 if(!trk) AliFatal("Not a standard AOD");
1824
d9117245 1825 if(trk->Pt() > 1) nHighPt++;
1826 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1827 qLepton[nLepton] = trk->Charge();
1828 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1829 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1830 mass[nLepton] = 0;
f052ef6f 1831 }
d9117245 1832 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1833 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1834 mass[nLepton] = 1;
1835 }
1836 nLepton++;
1837 }
1838 if(nLepton == 2){
1839 if(qLepton[0]*qLepton[1] < 0 && nHighPt > 0 && (mass[0]!=-1 || mass[1]!=-1)){
1840 vCandidate = vLepton[0]+vLepton[1];
1841 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1842 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1843 else {
1844 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1845 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
1846 }
1847 if(fChannel == -1 && vCandidate.Pt()<0.15) ((TH1D*)(fListJPsiTight->At(i)))->Fill(vCandidate.M());
1848 if(fChannel == 1 && vCandidate.Pt()<0.3) ((TH1D*)(fListJPsiTight->At(i)))->Fill(vCandidate.M());
f052ef6f 1849 }
1850 }
1851 }
d9117245 1852}//tight cuts
f052ef6f 1853
d9117245 1854//---------------------------------------------Psi2s------------------------------------------------------------------------
f052ef6f 1855
d9117245 1856 Double_t fPsi2sSels[4];
3d16cd00 1857
d9117245 1858 fPsi2sSels[0] = 50; //min number of TPC clusters
1859 fPsi2sSels[1] = 4; //chi2
1860 fPsi2sSels[2] = 2; //DCAz
1861 fPsi2sSels[3] = 4; // DCAxy 1x
3d16cd00 1862
d9117245 1863 Double_t fPsi2sSelsMid[4];
3d16cd00 1864
d9117245 1865 fPsi2sSelsMid[0] = 50; //min number of TPC clusters
1866 fPsi2sSelsMid[1] = 4; //chi2
1867 fPsi2sSelsMid[2] = 2; //DCAz
1868 fPsi2sSelsMid[3] = 4; // DCAxy 1x
03cff2a0 1869
d9117245 1870 Double_t fPsi2sSelsLoose[4];
3d16cd00 1871
d9117245 1872 fPsi2sSelsLoose[0] = 60; //min number of TPC clusters
1873 fPsi2sSelsLoose[1] = 5; //chi2
1874 fPsi2sSelsLoose[2] = 3; //DCAz
1875 fPsi2sSelsLoose[3] = 6; // DCAxy 2x
3d16cd00 1876
d9117245 1877 Double_t fPsi2sSelsTight[4];
3d16cd00 1878
d9117245 1879 fPsi2sSelsTight[0] = 70; //min number of TPC clusters
1880 fPsi2sSelsTight[1] = 3.5; //chi2
1881 fPsi2sSelsTight[2] = 1; //DCAz
1882 fPsi2sSelsTight[3] = 2; // DCAxy 0.5x
489951cd 1883
d9117245 1884 nGoodTracks = 0; nLepton=0; nHighPt=0; fChannel = 0;
1885 Int_t nSpdHits = 0;
1886 Double_t TrackPt[5]={0,0,0,0,0};
1887 Double_t MeanPt = -1;
489951cd 1888
d9117245 1889for(Int_t i=0; i<5; i++){
1890 //cout<<"Loose systematics psi2s, cut"<<i<<endl;
1891 for(Int_t j=0; j<4; j++){
1892 if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
1893 else fJPsiSels[j] = fJPsiSelsMid[j];
1894 }
1895
1896 //Four track loop
1897 nGoodTracks = 0; nSpdHits = 0;
1898 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 1899 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
3d16cd00 1900 if( !trk ) continue;
d9117245 1901 if(!(trk->TestFilterBit(1<<0))) continue;
3d16cd00 1902
1903 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1904 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
d9117245 1905 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1906 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1907 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1908 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1909 delete trk_clone;
1910 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
3d16cd00 1911
d9117245 1912 if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1913 if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1914 if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1915 if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1916 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1917
3d16cd00 1918 TrackIndex[nGoodTracks] = itr;
d9117245 1919 TrackPt[nGoodTracks] = trk->Pt();
3d16cd00 1920 nGoodTracks++;
d9117245 1921
1922 if(nGoodTracks > 4) break;
3d16cd00 1923 }//Track loop
d9117245 1924
1925 Int_t mass[3]={-1,-1,-1};
1926 fChannel = 0;
1927 nLepton=0; nPion=0; nHighPt=0;
3d16cd00 1928
d9117245 1929 if(nGoodTracks == 4){
1930 if(i!=4){ if(nSpdHits<2) continue;}
1931 MeanPt = GetMedian(TrackPt);
1932 for(Int_t k=0; k<4; k++){
f15c1f69 1933 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[k]));
1934 if(!trk) AliFatal("Not a standard AOD");
1935
d9117245 1936 if(trk->Pt() > MeanPt){
1937 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1938 qLepton[nLepton] = trk->Charge();
1939 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1940 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1941 mass[nLepton] = 0;
1942 }
1943 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1944 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1945 mass[nLepton] = 1;
1946 }
1947 nLepton++;
1948 }
1949 else{
1950 qPion[nPion] = trk->Charge();
1951 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1952 nPion++;
1953 }
1954 }
1955 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0) && mass[0] != -1 && mass[1] != -1){
1956 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1957 vDilepton = vLepton[0]+vLepton[1];
1958 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1959 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1960 else {
1961 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1962 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
1963 }
1964 if(fChannel == -1) if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) ((TH1D*)(fListPsi2sLoose->At(i)))->Fill(vCandidate.M());
1965 if(fChannel == 1) if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) ((TH1D*)(fListPsi2sLoose->At(i)))->Fill(vCandidate.M());
1966 }
1967 }
1968}//loose cuts
1969
1970for(Int_t i=0; i<4; i++){
1971 //cout<<"Tight systematics psi2s, cut"<<i<<endl;
1972 for(Int_t j=0; j<4; j++){
1973 if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
1974 else fJPsiSels[j] = fJPsiSelsMid[j];
1975 }
1976
03cff2a0 1977 //Four track loop
d9117245 1978 nGoodTracks = 0; nSpdHits = 0;
1979 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
f15c1f69 1980 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
03cff2a0 1981 if( !trk ) continue;
d9117245 1982 if(!(trk->TestFilterBit(1<<0))) continue;
03cff2a0
MB		 1983
1984 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1985 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
d9117245 1986 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1987 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1988 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1989 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1990 delete trk_clone;
1991 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
03cff2a0 1992
d9117245 1993 if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1994 if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1995 if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1996 if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1997 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1998
03cff2a0 1999 TrackIndex[nGoodTracks] = itr;
d9117245 2000 TrackPt[nGoodTracks] = trk->Pt();
03cff2a0 2001 nGoodTracks++;
d9117245 2002
2003 if(nGoodTracks > 4) break;
03cff2a0 2004 }//Track loop
d9117245 2005
2006 Int_t mass[3]={-1,-1,-1};
2007 fChannel = 0;
2008 nLepton=0; nPion=0; nHighPt=0;
03cff2a0 2009
3d16cd00 2010 if(nGoodTracks == 4){
d9117245 2011 if(nSpdHits<2) continue;
2012 MeanPt = GetMedian(TrackPt);
2013 for(Int_t k=0; k<4; k++){
f15c1f69 2014 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[k]));
2015 if(!trk) AliFatal("Not a standard AOD");
2016
d9117245 2017 if(trk->Pt() > MeanPt){
2018 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
2019 qLepton[nLepton] = trk->Charge();
2020 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
2021 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
2022 mass[nLepton] = 0;
2023 }
2024 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
2025 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
2026 mass[nLepton] = 1;
2027 }
2028 nLepton++;
2029 }
2030 else{
2031 qPion[nPion] = trk->Charge();
2032 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
2033 nPion++;
2034 }
2035 }
2036 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0) && mass[0] != -1 && mass[1] != -1){
2037 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
2038 vDilepton = vLepton[0]+vLepton[1];
2039 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
2040 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
2041 else {
2042 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
2043 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
2044 }
2045 if(fChannel == -1) if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) ((TH1D*)(fListPsi2sTight->At(i)))->Fill(vCandidate.M());
2046 if(fChannel == 1) if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) ((TH1D*)(fListPsi2sTight->At(i)))->Fill(vCandidate.M());
2047 }
2048 }
2049}//Tight cuts
22b569e5 2050
22b569e5 2051}