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