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