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