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