1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
23 #include "TClonesArray.h"
24 #include "TParticle.h"
25 #include "TObjString.h"
27 #include "TDatabasePDG.h"
28 #include "TLorentzVector.h"
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"
50 #include "AliCentrality.h"
51 #include "AliKFVertex.h"
52 #include "AliExternalTrackParam.h"
53 #include "AliTriggerAnalysis.h"
56 #include "AliAnalysisTaskUpcPsi2s.h"
58 ClassImp(AliAnalysisTaskUpcPsi2s);
63 //trees for UPC analysis,
64 // michal.broz@cern.ch
66 //_____________________________________________________________________________
67 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s()
68 : AliAnalysisTaskSE(),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fRunSystematics(kFALSE),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
69 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
70 fTOFtrig1(0), fTOFtrig2(0),
71 fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
72 fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
73 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
74 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
75 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
76 fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),fHistCint1TriggersPerRun(0),
77 fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
78 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),fHistDiLeptonMass(0),
79 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
80 fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0),fListPsi2sLoose(0),fListPsi2sTight(0)
86 }//AliAnalysisTaskUpcPsi2s
89 //_____________________________________________________________________________
90 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
91 : AliAnalysisTaskSE(name),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fRunSystematics(kFALSE),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
92 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
93 fTOFtrig1(0), fTOFtrig2(0),
94 fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
95 fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
96 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
97 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
98 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
99 fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),fHistCint1TriggersPerRun(0),
100 fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
101 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),fHistDiLeptonMass(0),
102 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
103 fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0),fListPsi2sLoose(0),fListPsi2sTight(0)
108 if( strstr(name,"ESD") ) fType = 0;
109 if( strstr(name,"AOD") ) fType = 1;
113 DefineOutput(1, TTree::Class());
114 DefineOutput(2, TTree::Class());
115 DefineOutput(3, TList::Class());
116 DefineOutput(4, TList::Class());
118 }//AliAnalysisTaskUpcPsi2s
120 //_____________________________________________________________________________
121 void AliAnalysisTaskUpcPsi2s::Init()
124 for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
125 for(Int_t i=0; i<4; i++) {
127 fPIDTPCMuon[i] = -666;
128 fPIDTPCElectron[i] = -666;
129 fPIDTPCPion[i] = -666;
130 fPIDTPCKaon[i] = -666;
131 fPIDTPCProton[i] = -666;
133 fPIDTOFMuon[i] = -666;
134 fPIDTOFElectron[i] = -666;
135 fPIDTOFPion[i] = -666;
136 fPIDTOFKaon[i] = -666;
137 fPIDTOFProton[i] = -666;
139 fTriggerInputsMC[i] = kFALSE;
141 for(Int_t i=0; i<3; i++){
145 fSpdVtxPos[i] = -666;
150 //_____________________________________________________________________________
151 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
171 }//~AliAnalysisTaskUpcPsi2s
174 //_____________________________________________________________________________
175 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
178 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
179 AliInputEventHandler *inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
180 fPIDResponse = inputHandler->GetPIDResponse();
183 fDataFilnam = new TObjString();
184 fDataFilnam->SetString("");
187 fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
188 fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
189 fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
190 fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
191 fGenPart = new TClonesArray("TParticle", 1000);
193 //output tree with JPsi candidate events
194 fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
195 fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
196 fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
197 fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
199 fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
200 fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
201 fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
202 fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
203 fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
204 fJPsiTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
205 fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
207 fJPsiTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
208 fJPsiTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
209 fJPsiTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
211 fJPsiTree ->Branch("fPIDTPCMuon", &fPIDTPCMuon[0], "fPIDTPCMuon[4]/D");
212 fJPsiTree ->Branch("fPIDTPCElectron", &fPIDTPCElectron[0], "fPIDTPCElectron[4]/D");
213 fJPsiTree ->Branch("fPIDTPCPion", &fPIDTPCPion[0], "fPIDTPCPion[4]/D");
214 fJPsiTree ->Branch("fPIDTPCKaon", &fPIDTPCKaon[0], "fPIDTPCKaon[4]/D");
215 fJPsiTree ->Branch("fPIDTPCProton", &fPIDTPCProton[0], "fPIDTPCProton[4]/D");
217 fJPsiTree ->Branch("fPIDTOFMuon", &fPIDTOFMuon[0], "fPIDTOFMuon[4]/D");
218 fJPsiTree ->Branch("fPIDTOFElectron", &fPIDTOFElectron[0], "fPIDTOFElectron[4]/D");
219 fJPsiTree ->Branch("fPIDTOFPion", &fPIDTOFPion[0], "fPIDTOFPion[4]/D");
220 fJPsiTree ->Branch("fPIDTOFKaon", &fPIDTOFKaon[0], "fPIDTOFKaon[4]/D");
221 fJPsiTree ->Branch("fPIDTOFProton", &fPIDTOFProton[0], "fPIDTOFProton[4]/D");
223 fJPsiTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
224 fJPsiTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
225 fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
226 fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
228 fJPsiTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
229 fJPsiTree ->Branch("fSpdVtxPos", &fSpdVtxPos[0], "fSpdVtxPos[3]/D");
231 fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
232 fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
233 fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
234 fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
235 fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
236 fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
237 fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
239 fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
242 fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
245 fJPsiTree ->Branch("fGenPart", &fGenPart);
246 fJPsiTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
250 //output tree with Psi2s candidate events
251 fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
252 fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
253 fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
254 fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
256 fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
257 fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
258 fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
259 fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
260 fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
261 fPsi2sTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
262 fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
264 fPsi2sTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
265 fPsi2sTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
266 fPsi2sTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
268 fPsi2sTree ->Branch("fPIDTPCMuon", &fPIDTPCMuon[0], "fPIDTPCMuon[4]/D");
269 fPsi2sTree ->Branch("fPIDTPCElectron", &fPIDTPCElectron[0], "fPIDTPCElectron[4]/D");
270 fPsi2sTree ->Branch("fPIDTPCPion", &fPIDTPCPion[0], "fPIDTPCPion[4]/D");
271 fPsi2sTree ->Branch("fPIDTPCKaon", &fPIDTPCKaon[0], "fPIDTPCKaon[4]/D");
272 fPsi2sTree ->Branch("fPIDTPCProton", &fPIDTPCProton[0], "fPIDTPCProton[4]/D");
274 fPsi2sTree ->Branch("fPIDTOFMuon", &fPIDTOFMuon[0], "fPIDTOFMuon[4]/D");
275 fPsi2sTree ->Branch("fPIDTOFElectron", &fPIDTOFElectron[0], "fPIDTOFElectron[4]/D");
276 fPsi2sTree ->Branch("fPIDTOFPion", &fPIDTOFPion[0], "fPIDTOFPion[4]/D");
277 fPsi2sTree ->Branch("fPIDTOFKaon", &fPIDTOFKaon[0], "fPIDTOFKaon[4]/D");
278 fPsi2sTree ->Branch("fPIDTOFProton", &fPIDTOFProton[0], "fPIDTOFProton[4]/D");
280 fPsi2sTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
281 fPsi2sTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
282 fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
283 fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
285 fPsi2sTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
286 fPsi2sTree ->Branch("fSpdVtxPos", &fSpdVtxPos[0], "fSpdVtxPos[3]/D");
288 fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
289 fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
290 fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
291 fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
292 fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
293 fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
294 fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
296 fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
299 fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
302 fPsi2sTree ->Branch("fGenPart", &fGenPart);
303 fPsi2sTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
307 fListTrig = new TList();
308 fListTrig ->SetOwner();
310 fHistCcup4TriggersPerRun = new TH1D("fHistCcup4TriggersPerRun", "fHistCcup4TriggersPerRun", 33000, 167000.5, 200000.5);
311 fListTrig->Add(fHistCcup4TriggersPerRun);
313 fHistCcup7TriggersPerRun = new TH1D("fHistCcup7TriggersPerRun", "fHistCcup7TriggersPerRun", 33000, 167000.5, 200000.5);
314 fListTrig->Add(fHistCcup7TriggersPerRun);
316 fHistCcup2TriggersPerRun = new TH1D("fHistCcup2TriggersPerRun", "fHistCcup2TriggersPerRun", 33000, 167000.5, 200000.5);
317 fListTrig->Add(fHistCcup2TriggersPerRun);
319 fHistCint1TriggersPerRun = new TH1D("fHistCint1TriggersPerRun", "fHistCint1TriggersPerRun", 33000, 167000.5, 200000.5);
320 fListTrig->Add(fHistCint1TriggersPerRun);
322 fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 33000, 167000.5, 200000.5);
323 fListTrig->Add(fHistZedTriggersPerRun);
325 fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 33000, 167000.5, 200000.5);
326 fListTrig->Add(fHistCvlnTriggersPerRun);
328 fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 33000, 167000.5, 200000.5);
329 fListTrig->Add(fHistMBTriggersPerRun);
331 fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 33000, 167000.5, 200000.5);
332 fListTrig->Add(fHistCentralTriggersPerRun);
334 fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 33000, 167000.5, 200000.5);
335 fListTrig->Add(fHistSemiCentralTriggersPerRun);
337 fListHist = new TList();
338 fListHist ->SetOwner();
340 TString CutNameJPsi[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Two good tracks",
341 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
342 fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",13,0.5,13.5);
343 for (Int_t i = 0; i<13; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
344 fListHist->Add(fHistNeventsJPsi);
346 fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
347 fListHist->Add(fHistTPCsignalJPsi);
349 fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
350 fListHist->Add(fHistDiLeptonPtJPsi);
352 fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
353 fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
354 fListHist->Add(fHistDiElectronMass);
356 fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
357 fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
358 fListHist->Add(fHistDiMuonMass);
360 fHistDiLeptonMass = new TH1D("fHistDiLeptonMass","Invariant mass of J/#psi candidates",130,2.1,6.0);
361 fHistDiLeptonMass->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}) (GeV/c)");
362 fListHist->Add(fHistDiLeptonMass);
364 TString CutNamePsi2s[14] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Four good tracks",
365 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
367 fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",14,0.5,14.5);
368 for (Int_t i = 0; i<14; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
369 fListHist->Add(fHistNeventsPsi2s);
371 fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
372 fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
373 fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
374 fListHist->Add(fHistPsi2sMassVsPt);
376 fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",50,2.5,5.5);
377 fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
378 fListHist->Add(fHistPsi2sMassCoherent);
380 fListSystematics = new TList();
381 fListSystematics->SetOwner();
382 fListSystematics->SetName("fListSystematics");
383 fListHist->Add(fListSystematics);
387 PostData(1, fJPsiTree);
388 PostData(2, fPsi2sTree);
389 PostData(3, fListTrig);
390 PostData(4, fListHist);
392 }//UserCreateOutputObjects
394 //_____________________________________________________________________________
395 void AliAnalysisTaskUpcPsi2s::InitSystematics()
398 fListJPsiLoose = new TList();
399 fListJPsiLoose->SetOwner();
400 fListJPsiLoose->SetName("JPsiLoose");
401 fListSystematics->Add(fListJPsiLoose);
403 TH1D *fHistJPsiNClusLoose = new TH1D("JPsiNClusLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
404 fListJPsiLoose->Add(fHistJPsiNClusLoose);
406 TH1D *fHistJPsiChi2Loose = new TH1D("JPsiChi2Loose","Invariant mass of J/#psi candidates",130,2.1,6.0);
407 fListJPsiLoose->Add(fHistJPsiChi2Loose);
409 TH1D *fHistJPsiDCAzLoose = new TH1D("JPsiDCAzLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
410 fListJPsiLoose->Add(fHistJPsiDCAzLoose);
412 TH1D *fHistJPsiDCAxyLoose = new TH1D("JPsiDCAxyLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
413 fListJPsiLoose->Add(fHistJPsiDCAxyLoose);
415 TH1D *fHistJPsiITShitsLoose = new TH1D("JPsiITShitsLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
416 fListJPsiLoose->Add(fHistJPsiITShitsLoose);
419 fListJPsiTight = new TList();
420 fListJPsiTight->SetOwner();
421 fListJPsiTight->SetName("JPsiTight");
422 fListSystematics->Add(fListJPsiTight);
424 TH1D *fHistJPsiNClusTight = new TH1D("JPsiNClusTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
425 fListJPsiTight->Add(fHistJPsiNClusTight);
427 TH1D *fHistJPsiChi2Tight = new TH1D("JPsiChi2Tight","Invariant mass of J/#psi candidates",130,2.1,6.0);
428 fListJPsiTight->Add(fHistJPsiChi2Tight);
430 TH1D *fHistJPsiDCAzTight = new TH1D("JPsiDCAzTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
431 fListJPsiTight->Add(fHistJPsiDCAzTight);
433 TH1D *fHistJPsiDCAxyTight = new TH1D("JPsiDCAxyTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
434 fListJPsiTight->Add(fHistJPsiDCAxyTight);
437 fListPsi2sLoose = new TList();
438 fListPsi2sLoose->SetOwner();
439 fListPsi2sLoose->SetName("Psi2sLoose");
440 fListSystematics->Add(fListPsi2sLoose);
442 TH1D *fHistPsi2sNClusLoose = new TH1D("Psi2sNClusLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
443 fListPsi2sLoose->Add(fHistPsi2sNClusLoose);
445 TH1D *fHistPsi2sChi2Loose = new TH1D("Psi2sChi2Loose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
446 fListPsi2sLoose->Add(fHistPsi2sChi2Loose);
448 TH1D *fHistPsi2sDCAzLoose = new TH1D("Psi2sDCAzLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
449 fListPsi2sLoose->Add(fHistPsi2sDCAzLoose);
451 TH1D *fHistPsi2sDCAxyLoose = new TH1D("Psi2sDCAxyLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
452 fListPsi2sLoose->Add(fHistPsi2sDCAxyLoose);
454 TH1D *fHistPsi2sITShitsLoose = new TH1D("Psi2sITShitsLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
455 fListPsi2sLoose->Add(fHistPsi2sITShitsLoose);
458 fListPsi2sTight = new TList();
459 fListPsi2sTight->SetOwner();
460 fListPsi2sTight->SetName("Psi2sTight");
461 fListSystematics->Add(fListPsi2sTight);
463 TH1D *fHistPsi2sNClusTight = new TH1D("Psi2sNClusTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
464 fListPsi2sTight->Add(fHistPsi2sNClusTight);
466 TH1D *fHistPsi2sChi2Tight = new TH1D("Psi2sChi2Tight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
467 fListPsi2sTight->Add(fHistPsi2sChi2Tight);
469 TH1D *fHistPsi2sDCAzTight = new TH1D("Psi2sDCAzTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
470 fListPsi2sTight->Add(fHistPsi2sDCAzTight);
472 TH1D *fHistPsi2sDCAxyTight = new TH1D("Psi2sDCAxyTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
473 fListPsi2sTight->Add(fHistPsi2sDCAxyTight);
478 //_____________________________________________________________________________
479 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
482 //cout<<"#################### Next event ##################"<<endl;
486 if(fRunHist) RunESDhist();
487 if(fRunTree) RunESDtree();
492 if(fRunHist) RunAODhist();
493 if(fRunTree) RunAODtree();
497 //_____________________________________________________________________________
498 void AliAnalysisTaskUpcPsi2s::RunAODtrig()
502 AliAODEvent *aod = (AliAODEvent*) InputEvent();
505 fRunNum = aod ->GetRunNumber();
507 TString trigger = aod->GetFiredTriggerClasses();
509 if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
510 if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
511 if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
513 if(trigger.Contains("CINT1")) fHistCint1TriggersPerRun->Fill(fRunNum); //CINT1 triggers
515 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
516 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
518 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
519 if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
521 //MB, Central and SemiCentral triggers
522 AliCentrality *centrality = aod->GetCentrality();
523 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
525 Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
526 //Double_t percentile = centrality->GetCentralityPercentile("V0M");
528 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
530 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
532 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
534 PostData(3, fListTrig);
537 //_____________________________________________________________________________
538 void AliAnalysisTaskUpcPsi2s::RunAODhist()
541 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
543 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
544 Double_t muonMass = partMuon->Mass();
546 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
547 Double_t electronMass = partElectron->Mass();
549 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
550 Double_t pionMass = partPion->Mass();
553 AliAODEvent *aod = (AliAODEvent*) InputEvent();
556 // cout<<"Event number: "<<((TTree*) GetInputData(0))->GetTree()->GetReadEntry()<<endl;
558 fHistNeventsJPsi->Fill(1);
559 fHistNeventsPsi2s->Fill(1);
562 TString trigger = aod->GetFiredTriggerClasses();
564 if(!isMC && !trigger.Contains("CCUP") ) return;
566 fHistNeventsJPsi->Fill(2);
567 fHistNeventsPsi2s->Fill(2);
570 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
571 fVtxContrib = fAODVertex->GetNContributors();
572 if(fVtxContrib < 2) return;
574 fHistNeventsJPsi->Fill(3);
575 fHistNeventsPsi2s->Fill(3);
578 AliAODVZERO *fV0data = aod ->GetVZEROData();
579 AliAODZDC *fZDCdata = aod->GetZDCData();
581 fV0Adecision = fV0data->GetV0ADecision();
582 fV0Cdecision = fV0data->GetV0CDecision();
583 if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
585 fHistNeventsJPsi->Fill(4);
586 fHistNeventsPsi2s->Fill(4);
588 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
589 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
591 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
593 fHistNeventsJPsi->Fill(5);
594 fHistNeventsPsi2s->Fill(5);
596 //Systematics - cut variation
597 if(fRunSystematics) RunAODsystematics(aod);
600 Int_t nGoodTracks = 0;
601 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
603 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
604 Short_t qLepton[4], qPion[4];
605 UInt_t nLepton=0, nPion=0, nHighPt=0, nSpdHits=0;
606 Double_t fRecTPCsignal[5], fRecTPCsignalDist;
608 Int_t mass[3]={-1,-1,-1};
609 Double_t TrackPt[5]={0,0,0,0,0};
610 Double_t MeanPt = -1;
614 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
615 AliAODTrack *trk = aod->GetTrack(itr);
617 if(!(trk->TestFilterBit(1<<0))) continue;
619 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
620 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
621 if(trk->GetTPCNcls() < 50)continue;
622 if(trk->Chi2perNDF() > 4)continue;
623 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
624 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
625 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
627 if(TMath::Abs(dca[1]) > 2) continue;
628 Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
629 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
630 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
632 TrackIndex[nGoodTracks] = itr;
633 TrackPt[nGoodTracks] = trk->Pt();
636 if(nGoodTracks > 4) break;
639 nLepton=0; nPion=0; nHighPt=0;
640 mass[0]= -1; mass[1]= -1, mass[2]= -1;
642 if(nGoodTracks == 4 && nSpdHits>1){
643 MeanPt = GetMedian(TrackPt);
644 fHistNeventsPsi2s->Fill(6);
645 for(Int_t i=0; i<4; i++){
646 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
648 if(trk->Pt() > MeanPt){
649 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
650 qLepton[nLepton] = trk->Charge();
651 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
652 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
655 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
656 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
662 qPion[nPion] = trk->Charge();
663 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
667 if(nLepton > 2 || nPion > 2) break;
669 if((nLepton == 2) && (nPion == 2)){
670 fHistNeventsPsi2s->Fill(7);
671 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(8);
672 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(9);
673 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
674 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
675 fHistNeventsPsi2s->Fill(11);
676 if(mass[0] != -1 && mass[1] != -1) {
677 fHistNeventsPsi2s->Fill(12);
678 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
679 vDilepton = vLepton[0]+vLepton[1];
680 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
681 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
682 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
684 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
685 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
689 fHistNeventsPsi2s->Fill(13);
690 if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
693 fHistNeventsPsi2s->Fill(14);
694 if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) fHistPsi2sMassCoherent->Fill(vCandidate.M());
703 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
704 AliAODTrack *trk = aod->GetTrack(itr);
706 if(!(trk->TestFilterBit(1<<0))) continue;
708 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
709 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
710 if(trk->GetTPCNcls() < 70)continue;
711 if(trk->Chi2perNDF() > 4)continue;
712 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
713 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
714 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
715 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
717 if(TMath::Abs(dca[1]) > 2) continue;
718 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
719 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
721 TrackIndex[nGoodTracks] = itr;
724 if(nGoodTracks > 2) break;
727 nLepton=0; nPion=0; nHighPt=0;
728 mass[0]= -1; mass[1]= -1, mass[2]= -1;
730 if(nGoodTracks == 2){
731 fHistNeventsJPsi->Fill(6);
732 for(Int_t i=0; i<2; i++){
733 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
734 if(trk->Pt() > 1) nHighPt++;
735 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
736 qLepton[nLepton] = trk->Charge();
737 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
738 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
741 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
742 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
748 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(7);
749 if(qLepton[0]*qLepton[1] < 0){
750 fHistNeventsJPsi->Fill(8);
752 fHistNeventsJPsi->Fill(9);
753 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
754 if(nHighPt == 2) fHistNeventsJPsi->Fill(10);
755 if(mass[0] != -1 && mass[1] != -1) {
756 fHistNeventsJPsi->Fill(11);
757 vCandidate = vLepton[0]+vLepton[1];
758 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
759 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
761 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
762 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
764 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
766 fHistDiMuonMass->Fill(vCandidate.M());
767 if(vCandidate.Pt()<0.15)fHistDiLeptonMass->Fill(vCandidate.M());
768 fHistNeventsJPsi->Fill(12);
771 fHistDiElectronMass->Fill(vCandidate.M());
772 if(vCandidate.Pt()<0.3)fHistDiLeptonMass->Fill(vCandidate.M());
773 fHistNeventsJPsi->Fill(13);
782 PostData(4, fListHist);
786 //_____________________________________________________________________________
787 void AliAnalysisTaskUpcPsi2s::RunAODtree()
790 AliAODEvent *aod = (AliAODEvent*) InputEvent();
793 if(isMC) RunAODMC(aod);
796 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
797 fDataFilnam->Clear();
798 fDataFilnam->SetString(filnam);
799 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
800 fRunNum = aod ->GetRunNumber();
803 TString trigger = aod->GetFiredTriggerClasses();
805 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
806 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
807 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
808 fTrigger[3] = trigger.Contains("CINT1");
810 Bool_t isTriggered = kFALSE;
811 for(Int_t i=0; i<ntrg; i++) {
812 if( fTrigger[i] ) isTriggered = kTRUE;
814 if(!isMC && !isTriggered ) return;
817 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
818 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
820 //Event identification
821 fPerNum = aod ->GetPeriodNumber();
822 fOrbNum = aod ->GetOrbitNumber();
823 fBCrossNum = aod ->GetBunchCrossNumber();
826 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
827 fVtxContrib = fAODVertex->GetNContributors();
828 fVtxPos[0] = fAODVertex->GetX();
829 fVtxPos[1] = fAODVertex->GetY();
830 fVtxPos[2] = fAODVertex->GetZ();
832 fAODVertex->GetCovarianceMatrix(CovMatx);
833 fVtxErr[0] = CovMatx[0];
834 fVtxErr[1] = CovMatx[1];
835 fVtxErr[2] = CovMatx[2];
836 fVtxChi2 = fAODVertex->GetChi2();
837 fVtxNDF = fAODVertex->GetNDF();
840 AliAODVertex *fSPDVertex = aod->GetPrimaryVertexSPD();
842 fSpdVtxPos[0] = fSPDVertex->GetX();
843 fSpdVtxPos[1] = fSPDVertex->GetY();
844 fSpdVtxPos[2] = fSPDVertex->GetZ();
848 fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
851 AliAODVZERO *fV0data = aod ->GetVZEROData();
852 AliAODZDC *fZDCdata = aod->GetZDCData();
854 fV0Adecision = fV0data->GetV0ADecision();
855 fV0Cdecision = fV0data->GetV0CDecision();
856 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
857 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
861 //Track loop - loose cuts
862 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
863 AliAODTrack *trk = aod->GetTrack(itr);
865 if(!(trk->TestFilterBit(1<<0))) continue;
867 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
868 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
869 if(trk->GetTPCNcls() < 20)continue;
871 }//Track loop -loose cuts
874 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
877 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
878 AliAODTrack *trk = aod->GetTrack(itr);
880 if(!(trk->TestFilterBit(1<<0))) continue;
882 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
883 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
884 if(trk->GetTPCNcls() < 70)continue;
885 if(trk->Chi2perNDF() > 4)continue;
886 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
887 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
888 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
889 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
891 if(TMath::Abs(dca[1]) > 2) continue;
892 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
893 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
896 TrackIndex[nGoodTracks] = itr;
899 if(nGoodTracks > 2) break;
902 fJPsiAODTracks->Clear("C");
903 if(nGoodTracks == 2){
905 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
906 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
907 Double_t muonMass = partMuon->Mass();
908 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
909 Double_t electronMass = partElectron->Mass();
910 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
911 Double_t pionMass = partPion->Mass();
915 Double_t KFmass = pionMass;
916 Double_t fRecTPCsignal;
917 AliKFParticle *KFpart[2];
918 AliKFVertex *KFvtx = new AliKFVertex();
919 KFvtx->SetField(aod->GetMagneticField());
921 for(Int_t i=0; i<2; i++){
922 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
924 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
925 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
926 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
929 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
930 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
932 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
933 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
934 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
935 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
936 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
938 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
939 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
940 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
941 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
942 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
944 trk->GetPosition(KFpar);
945 trk->PxPyPz(KFpar+3);
946 trk->GetCovMatrix(KFcov);
949 fRecTPCsignal = trk->GetTPCsignal();
950 if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
951 if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
953 else KFmass = pionMass;
955 KFpart[i] = new AliKFParticle();
956 KFpart[i]->SetField(aod->GetMagneticField());
957 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
958 KFvtx->AddDaughter(*KFpart[i]);
961 Double_t pos[3]={0,0,0};
962 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
963 parTrk->CopyFromVTrack((AliVTrack*) trk);
964 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
966 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
967 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
971 fKfVtxPos[0]= KFvtx->GetX();
972 fKfVtxPos[1]= KFvtx->GetY();
973 fKfVtxPos[2]= KFvtx->GetZ();
974 for(UInt_t i=0; i<2; i++)delete KFpart[i];
977 if(!isMC) fJPsiTree ->Fill();
982 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
983 AliAODTrack *trk = aod->GetTrack(itr);
985 if(!(trk->TestFilterBit(1<<0))) continue;
987 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
988 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
989 if(trk->GetTPCNcls() < 50)continue;
990 if(trk->Chi2perNDF() > 4)continue;
991 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
992 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
993 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
995 if(TMath::Abs(dca[1]) > 2) continue;
996 Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
997 if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
999 TrackIndex[nGoodTracks] = itr;
1002 if(nGoodTracks > 4) break;
1005 fPsi2sAODTracks->Clear("C");
1006 if(nGoodTracks == 4){
1008 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1009 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1010 Double_t muonMass = partMuon->Mass();
1011 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1012 Double_t electronMass = partElectron->Mass();
1013 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1014 Double_t pionMass = partPion->Mass();
1018 Double_t KFmass = pionMass;
1019 Double_t fRecTPCsignal;
1020 AliKFParticle *KFpart[4];
1021 AliKFVertex *KFvtx = new AliKFVertex();
1022 KFvtx->SetField(aod->GetMagneticField());
1024 for(Int_t i=0; i<4; i++){
1025 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
1027 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1028 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1029 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1032 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
1033 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
1036 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1037 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1038 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1039 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1040 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
1042 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1043 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1044 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1045 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1046 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
1048 trk->GetPosition(KFpar);
1049 trk->PxPyPz(KFpar+3);
1050 trk->GetCovMatrix(KFcov);
1053 fRecTPCsignal = trk->GetTPCsignal();
1054 if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
1055 if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
1057 else KFmass = pionMass;
1059 KFpart[i] = new AliKFParticle();
1060 KFpart[i]->SetField(aod->GetMagneticField());
1061 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
1062 KFvtx->AddDaughter(*KFpart[i]);
1064 Double_t pos[3]={0,0,0};
1065 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
1066 parTrk->CopyFromVTrack((AliVTrack*) trk);
1067 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
1069 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1070 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1074 fKfVtxPos[0]= KFvtx->GetX();
1075 fKfVtxPos[1]= KFvtx->GetY();
1076 fKfVtxPos[2]= KFvtx->GetZ();
1077 for(UInt_t i=0; i<4; i++)delete KFpart[i];
1079 if(!isMC) fPsi2sTree ->Fill();
1084 fPsi2sTree ->Fill();
1087 PostData(1, fJPsiTree);
1088 PostData(2, fPsi2sTree);
1093 //_____________________________________________________________________________
1094 void AliAnalysisTaskUpcPsi2s::RunAODMC(AliAODEvent *aod)
1097 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
1098 fTriggerInputsMC[0] = kFALSE;//0SM2
1099 fTriggerInputsMC[1] = fL0inputs & (1 << 0);//0VBA
1100 fTriggerInputsMC[2] = fL0inputs & (1 << 1);//0VBC
1101 fTriggerInputsMC[3] = fL0inputs & (1 << 9);//0OMU
1103 fGenPart->Clear("C");
1105 TClonesArray *arrayMC = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1106 if(!arrayMC) return;
1109 //loop over mc particles
1110 for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
1111 AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(imc);
1112 if(!mcPart) continue;
1114 if(mcPart->GetMother() >= 0) continue;
1116 TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1117 part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1118 part->SetPdgCode(mcPart->GetPdgCode());
1119 part->SetUniqueID(imc);
1120 }//loop over mc particles
1125 //_____________________________________________________________________________
1126 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
1130 AliESDEvent *esd = (AliESDEvent*) InputEvent();
1133 fRunNum = esd ->GetRunNumber();
1135 TString trigger = esd->GetFiredTriggerClasses();
1137 if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
1138 if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
1139 if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
1141 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
1142 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
1144 if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
1146 //MB, Central and SemiCentral triggers
1147 AliCentrality *centrality = esd->GetCentrality();
1148 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1150 //Double_t percentile = centrality->GetCentralityPercentile("V0M");
1151 Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
1153 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
1155 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
1157 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
1160 PostData(3, fListTrig);
1163 //_____________________________________________________________________________
1164 void AliAnalysisTaskUpcPsi2s::RunESDhist()
1167 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1169 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1170 Double_t muonMass = partMuon->Mass();
1172 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1173 Double_t electronMass = partElectron->Mass();
1175 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1176 Double_t pionMass = partPion->Mass();
1179 AliESDEvent *esd = (AliESDEvent*) InputEvent();
1182 fHistNeventsJPsi->Fill(1);
1183 fHistNeventsPsi2s->Fill(1);
1186 TString trigger = esd->GetFiredTriggerClasses();
1188 if(!isMC && !trigger.Contains("CCUP") ) return;
1190 fHistNeventsJPsi->Fill(2);
1191 fHistNeventsPsi2s->Fill(2);
1194 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1195 fVtxContrib = fESDVertex->GetNContributors();
1196 if(fVtxContrib < 2) return;
1198 fHistNeventsJPsi->Fill(3);
1199 fHistNeventsPsi2s->Fill(3);
1202 AliESDVZERO *fV0data = esd->GetVZEROData();
1203 AliESDZDC *fZDCdata = esd->GetESDZDC();
1205 fV0Adecision = fV0data->GetV0ADecision();
1206 fV0Cdecision = fV0data->GetV0CDecision();
1207 if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
1209 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1210 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1211 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
1213 fHistNeventsJPsi->Fill(4);
1214 fHistNeventsPsi2s->Fill(4);
1216 Int_t nGoodTracks=0;
1218 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1220 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1221 Short_t qLepton[4], qPion[4];
1222 UInt_t nLepton=0, nPion=0, nHighPt=0;
1223 Double_t fRecTPCsignal[5];
1224 Int_t mass[3]={-1,-1,-1};
1227 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1228 AliESDtrack *trk = esd->GetTrack(itr);
1229 if( !trk ) continue;
1231 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1232 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1233 if(trk->GetTPCNcls() < 70)continue;
1234 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1235 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1236 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1237 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1238 trk->GetImpactParameters(dca[0],dca[1]);
1239 if(TMath::Abs(dca[1]) > 2) continue;
1240 if(TMath::Abs(dca[1]) > 0.2) continue;
1242 TrackIndex[nGoodTracks] = itr;
1244 if(nGoodTracks > 2) break;
1248 if(nGoodTracks == 2){
1249 fHistNeventsJPsi->Fill(5);
1250 for(Int_t i=0; i<2; i++){
1251 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1252 if(trk->Pt() > 1) nHighPt++;
1253 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1254 qLepton[nLepton] = trk->Charge();
1255 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1256 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1259 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1260 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1266 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
1267 if(qLepton[0]*qLepton[1] < 0){
1268 fHistNeventsJPsi->Fill(7);
1270 fHistNeventsJPsi->Fill(8);
1271 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
1272 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
1273 if(mass[0] == mass[1] && mass[0] != -1) {
1274 fHistNeventsJPsi->Fill(10);
1275 vCandidate = vLepton[0]+vLepton[1];
1276 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
1278 fHistDiMuonMass->Fill(vCandidate.M());
1279 fHistNeventsJPsi->Fill(11);
1282 fHistDiElectronMass->Fill(vCandidate.M());
1283 fHistNeventsJPsi->Fill(12);
1290 nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1291 mass[0]= -1; mass[1]= -1, mass[2]= -1;
1294 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1295 AliESDtrack *trk = esd->GetTrack(itr);
1296 if( !trk ) continue;
1298 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1299 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1300 if(trk->GetTPCNcls() < 50)continue;
1301 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1302 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1303 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1304 trk->GetImpactParameters(dca[0],dca[1]);
1305 if(TMath::Abs(dca[1]) > 2) continue;
1307 TrackIndex[nGoodTracks] = itr;
1309 if(nGoodTracks > 4) break;
1312 if(nGoodTracks == 4){
1313 fHistNeventsPsi2s->Fill(5);
1314 for(Int_t i=0; i<4; i++){
1315 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1318 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1319 qLepton[nLepton] = trk->Charge();
1320 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1321 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1324 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1325 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1331 qPion[nPion] = trk->Charge();
1332 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1336 if(nLepton > 2 || nPion > 2) break;
1338 if((nLepton == 2) && (nPion == 2)){
1339 fHistNeventsPsi2s->Fill(6);
1340 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1341 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1342 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1343 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1344 fHistNeventsPsi2s->Fill(10);
1345 if(mass[0] == mass[1]) {
1346 fHistNeventsPsi2s->Fill(11);
1347 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1348 vDilepton = vLepton[0]+vLepton[1];
1349 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1350 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1351 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1352 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1358 PostData(4, fListHist);
1362 //_____________________________________________________________________________
1363 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1367 AliESDEvent *esd = (AliESDEvent*) InputEvent();
1370 if(isMC) RunESDMC(esd);
1373 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1374 fDataFilnam->Clear();
1375 fDataFilnam->SetString(filnam);
1376 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1377 fRunNum = esd->GetRunNumber();
1380 TString trigger = esd->GetFiredTriggerClasses();
1382 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1383 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
1384 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1386 Bool_t isTriggered = kFALSE;
1387 for(Int_t i=0; i<ntrg; i++) {
1388 if( fTrigger[i] ) isTriggered = kTRUE;
1390 if(!isMC && !isTriggered ) return;
1393 fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1394 fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1396 //TOF trigger info (0OMU)
1397 fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1398 fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1400 //Event identification
1401 fPerNum = esd->GetPeriodNumber();
1402 fOrbNum = esd->GetOrbitNumber();
1403 fBCrossNum = esd->GetBunchCrossNumber();
1406 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1407 fVtxContrib = fESDVertex->GetNContributors();
1408 fVtxPos[0] = fESDVertex->GetX();
1409 fVtxPos[1] = fESDVertex->GetY();
1410 fVtxPos[2] = fESDVertex->GetZ();
1411 Double_t CovMatx[6];
1412 fESDVertex->GetCovarianceMatrix(CovMatx);
1413 fVtxErr[0] = CovMatx[0];
1414 fVtxErr[1] = CovMatx[1];
1415 fVtxErr[2] = CovMatx[2];
1416 fVtxChi2 = fESDVertex->GetChi2();
1417 fVtxNDF = fESDVertex->GetNDF();
1419 //SPD primary vertex
1420 AliESDVertex *fSPDVertex = (AliESDVertex*) esd->GetPrimaryVertexSPD();
1422 fSpdVtxPos[0] = fSPDVertex->GetX();
1423 fSpdVtxPos[1] = fSPDVertex->GetY();
1424 fSpdVtxPos[2] = fSPDVertex->GetZ();
1428 fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1431 AliESDVZERO *fV0data = esd->GetVZEROData();
1432 AliESDZDC *fZDCdata = esd->GetESDZDC();
1434 fV0Adecision = fV0data->GetV0ADecision();
1435 fV0Cdecision = fV0data->GetV0CDecision();
1436 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1437 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1441 //Track loop - loose cuts
1442 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1443 AliESDtrack *trk = esd->GetTrack(itr);
1444 if( !trk ) continue;
1446 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1447 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1448 if(trk->GetTPCNcls() < 20)continue;
1450 }//Track loop -loose cuts
1452 Int_t nGoodTracks=0;
1453 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1456 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1457 AliESDtrack *trk = esd->GetTrack(itr);
1458 if( !trk ) continue;
1460 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1461 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1462 if(trk->GetTPCNcls() < 70)continue;
1463 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1464 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1465 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1466 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1467 trk->GetImpactParameters(dca[0],dca[1]);
1468 if(TMath::Abs(dca[1]) > 2) continue;
1469 if(TMath::Abs(dca[0]) > 0.2) continue;
1471 TrackIndex[nGoodTracks] = itr;
1473 if(nGoodTracks > 2) break;
1476 fJPsiESDTracks->Clear("C");
1477 if(nGoodTracks == 2){
1478 for(Int_t i=0; i<2; i++){
1479 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1481 AliExternalTrackParam cParam;
1482 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1484 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
1486 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1487 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1488 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1489 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1490 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
1492 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1493 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1494 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1495 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1496 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
1498 Double_t pos[3]={0,0,0};
1499 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1501 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1502 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1505 if(!isMC) fJPsiTree ->Fill();
1510 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1511 AliESDtrack *trk = esd->GetTrack(itr);
1512 if( !trk ) continue;
1514 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1515 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1516 if(trk->GetTPCNcls() < 50)continue;
1517 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1518 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1519 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1520 trk->GetImpactParameters(dca[0],dca[1]);
1521 if(TMath::Abs(dca[1]) > 2) continue;
1523 TrackIndex[nGoodTracks] = itr;
1525 if(nGoodTracks > 4) break;
1528 fPsi2sESDTracks->Clear("C");
1529 if(nGoodTracks == 4){
1530 for(Int_t i=0; i<4; i++){
1531 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1533 AliExternalTrackParam cParam;
1534 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1536 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1538 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1539 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1540 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1541 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1542 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
1544 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1545 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1546 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1547 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1548 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
1550 Double_t pos[3]={0,0,0};
1551 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1553 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1554 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1557 if(!isMC) fPsi2sTree ->Fill();
1562 fPsi2sTree ->Fill();
1565 PostData(1, fJPsiTree);
1566 PostData(2, fPsi2sTree);
1571 //_____________________________________________________________________________
1572 void AliAnalysisTaskUpcPsi2s::RunESDMC(AliESDEvent* esd)
1575 AliTriggerAnalysis *fTrigAna = new AliTriggerAnalysis();
1576 fTrigAna->SetAnalyzeMC(isMC);
1578 if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) > 1) fTriggerInputsMC[0] = kTRUE;
1579 if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) < 2) fTriggerInputsMC[0] = kFALSE;
1581 fTriggerInputsMC[1] = esd->GetHeader()->IsTriggerInputFired("0OMU");
1582 fTriggerInputsMC[2] = esd->GetHeader()->IsTriggerInputFired("0VBA");
1583 fTriggerInputsMC[3] = esd->GetHeader()->IsTriggerInputFired("0VBC");
1585 fGenPart->Clear("C");
1587 AliMCEvent *mc = MCEvent();
1591 //loop over mc particles
1592 for(Int_t imc=0; imc<mc->GetNumberOfTracks(); imc++) {
1593 AliMCParticle *mcPart = (AliMCParticle*) mc->GetTrack(imc);
1594 if(!mcPart) continue;
1596 if(mcPart->GetMother() >= 0) continue;
1598 TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1599 part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1600 part->SetPdgCode(mcPart->PdgCode());
1601 part->SetUniqueID(imc);
1602 }//loop over mc particles
1608 //_____________________________________________________________________________
1609 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
1612 cout<<"Analysis complete."<<endl;
1615 //_____________________________________________________________________________
1616 Double_t AliAnalysisTaskUpcPsi2s::GetMedian(Double_t *daArray) {
1617 // Allocate an array of the same size and sort it.
1618 Double_t dpSorted[4];
1619 for (Int_t i = 0; i < 4; ++i) {
1620 dpSorted[i] = daArray[i];
1622 for (Int_t i = 3; i > 0; --i) {
1623 for (Int_t j = 0; j < i; ++j) {
1624 if (dpSorted[j] > dpSorted[j+1]) {
1625 Double_t dTemp = dpSorted[j];
1626 dpSorted[j] = dpSorted[j+1];
1627 dpSorted[j+1] = dTemp;
1632 // Middle or average of middle values in the sorted array.
1633 Double_t dMedian = 0.0;
1634 dMedian = (dpSorted[2] + dpSorted[1])/2.0;
1639 //_____________________________________________________________________________
1640 void AliAnalysisTaskUpcPsi2s::RunAODsystematics(AliAODEvent* aod)
1643 Double_t fJPsiSels[4];
1645 fJPsiSels[0] = 70; //min number of TPC clusters
1646 fJPsiSels[1] = 4; //chi2
1647 fJPsiSels[2] = 2; //DCAz
1648 fJPsiSels[3] = 1; // DCAxy 1x
1650 Double_t fJPsiSelsMid[4];
1652 fJPsiSelsMid[0] = 70; //min number of TPC clusters
1653 fJPsiSelsMid[1] = 4; //chi2
1654 fJPsiSelsMid[2] = 2; //DCAz
1655 fJPsiSelsMid[3] = 1; // DCAxy 1x
1657 Double_t fJPsiSelsLoose[4];
1659 fJPsiSelsLoose[0] = 60; //min number of TPC clusters
1660 fJPsiSelsLoose[1] = 5; //chi2
1661 fJPsiSelsLoose[2] = 3; //DCAz
1662 fJPsiSelsLoose[3] = 2; // DCAxy 2x
1664 Double_t fJPsiSelsTight[4];
1666 fJPsiSelsTight[0] = 80; //min number of TPC clusters
1667 fJPsiSelsTight[1] = 3.5; //chi2
1668 fJPsiSelsTight[2] = 1; //DCAz
1669 fJPsiSelsTight[3] = 0.5; // DCAxy 0.5x
1671 Int_t nGoodTracks = 0;
1672 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1674 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1675 Short_t qLepton[4],qPion[4];
1676 UInt_t nLepton=0, nPion=0, nHighPt=0;
1677 Double_t fRecTPCsignal[5], fRecTPCsignalDist;
1680 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
1682 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1684 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1685 Double_t muonMass = partMuon->Mass();
1687 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1688 Double_t electronMass = partElectron->Mass();
1690 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1691 Double_t pionMass = partPion->Mass();
1694 for(Int_t i=0; i<5; i++){
1695 //cout<<"Loose sytematics, cut"<<i<<endl;
1696 for(Int_t j=0; j<4; j++){
1697 if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
1698 else fJPsiSels[j] = fJPsiSelsMid[j];
1702 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1703 AliAODTrack *trk = aod->GetTrack(itr);
1704 if( !trk ) continue;
1705 if(!(trk->TestFilterBit(1<<0))) continue;
1707 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1708 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1709 if(i!=4){ if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;}
1710 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1711 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1712 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1714 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1716 if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1717 if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1718 if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1719 if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1721 TrackIndex[nGoodTracks] = itr;
1724 if(nGoodTracks > 2) break;
1727 Int_t mass[3]={-1,-1,-1};
1729 nLepton=0; nHighPt=0;
1731 if(nGoodTracks == 2){
1732 for(Int_t k=0; k<2; k++){
1733 AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);
1734 if(trk->Pt() > 1) nHighPt++;
1735 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1736 qLepton[nLepton] = trk->Charge();
1737 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1738 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1741 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1742 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1748 if(qLepton[0]*qLepton[1] < 0 && nHighPt > 0 && (mass[0]!=-1 || mass[1]!=-1)){
1749 vCandidate = vLepton[0]+vLepton[1];
1750 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1751 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1753 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1754 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
1756 if(fChannel == -1 && vCandidate.Pt()<0.15) ((TH1D*)(fListJPsiLoose->At(i)))->Fill(vCandidate.M());
1757 if(fChannel == 1 && vCandidate.Pt()<0.3) ((TH1D*)(fListJPsiLoose->At(i)))->Fill(vCandidate.M());
1763 for(Int_t i=0; i<4; i++){
1764 //cout<<"Tight sytematics, cut"<<i<<endl;
1765 for(Int_t j=0; j<4; j++){
1766 if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
1767 else fJPsiSels[j] = fJPsiSelsMid[j];
1771 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1772 AliAODTrack *trk = aod->GetTrack(itr);
1773 if( !trk ) continue;
1774 if(!(trk->TestFilterBit(1<<0))) continue;
1776 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1777 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1778 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1779 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1780 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1781 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1783 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1785 if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1786 if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1787 if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1788 if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1790 TrackIndex[nGoodTracks] = itr;
1793 if(nGoodTracks > 2) break;
1796 Int_t mass[3]={-1,-1,-1};
1798 nLepton=0; nHighPt=0;
1800 if(nGoodTracks == 2){
1801 for(Int_t k=0; k<2; k++){
1802 AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);
1803 if(trk->Pt() > 1) nHighPt++;
1804 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1805 qLepton[nLepton] = trk->Charge();
1806 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1807 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1810 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1811 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1817 if(qLepton[0]*qLepton[1] < 0 && nHighPt > 0 && (mass[0]!=-1 || mass[1]!=-1)){
1818 vCandidate = vLepton[0]+vLepton[1];
1819 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1820 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1822 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1823 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
1825 if(fChannel == -1 && vCandidate.Pt()<0.15) ((TH1D*)(fListJPsiTight->At(i)))->Fill(vCandidate.M());
1826 if(fChannel == 1 && vCandidate.Pt()<0.3) ((TH1D*)(fListJPsiTight->At(i)))->Fill(vCandidate.M());
1832 //---------------------------------------------Psi2s------------------------------------------------------------------------
1834 Double_t fPsi2sSels[4];
1836 fPsi2sSels[0] = 50; //min number of TPC clusters
1837 fPsi2sSels[1] = 4; //chi2
1838 fPsi2sSels[2] = 2; //DCAz
1839 fPsi2sSels[3] = 4; // DCAxy 1x
1841 Double_t fPsi2sSelsMid[4];
1843 fPsi2sSelsMid[0] = 50; //min number of TPC clusters
1844 fPsi2sSelsMid[1] = 4; //chi2
1845 fPsi2sSelsMid[2] = 2; //DCAz
1846 fPsi2sSelsMid[3] = 4; // DCAxy 1x
1848 Double_t fPsi2sSelsLoose[4];
1850 fPsi2sSelsLoose[0] = 60; //min number of TPC clusters
1851 fPsi2sSelsLoose[1] = 5; //chi2
1852 fPsi2sSelsLoose[2] = 3; //DCAz
1853 fPsi2sSelsLoose[3] = 6; // DCAxy 2x
1855 Double_t fPsi2sSelsTight[4];
1857 fPsi2sSelsTight[0] = 70; //min number of TPC clusters
1858 fPsi2sSelsTight[1] = 3.5; //chi2
1859 fPsi2sSelsTight[2] = 1; //DCAz
1860 fPsi2sSelsTight[3] = 2; // DCAxy 0.5x
1862 nGoodTracks = 0; nLepton=0; nHighPt=0; fChannel = 0;
1864 Double_t TrackPt[5]={0,0,0,0,0};
1865 Double_t MeanPt = -1;
1867 for(Int_t i=0; i<5; i++){
1868 //cout<<"Loose systematics psi2s, cut"<<i<<endl;
1869 for(Int_t j=0; j<4; j++){
1870 if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
1871 else fJPsiSels[j] = fJPsiSelsMid[j];
1875 nGoodTracks = 0; nSpdHits = 0;
1876 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1877 AliAODTrack *trk = aod->GetTrack(itr);
1878 if( !trk ) continue;
1879 if(!(trk->TestFilterBit(1<<0))) continue;
1881 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1882 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1883 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1884 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1885 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1886 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1888 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1890 if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1891 if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1892 if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1893 if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1894 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1896 TrackIndex[nGoodTracks] = itr;
1897 TrackPt[nGoodTracks] = trk->Pt();
1900 if(nGoodTracks > 4) break;
1903 Int_t mass[3]={-1,-1,-1};
1905 nLepton=0; nPion=0; nHighPt=0;
1907 if(nGoodTracks == 4){
1908 if(i!=4){ if(nSpdHits<2) continue;}
1909 MeanPt = GetMedian(TrackPt);
1910 for(Int_t k=0; k<4; k++){
1911 AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);
1912 if(trk->Pt() > MeanPt){
1913 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1914 qLepton[nLepton] = trk->Charge();
1915 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1916 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1919 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1920 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1926 qPion[nPion] = trk->Charge();
1927 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1931 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0) && mass[0] != -1 && mass[1] != -1){
1932 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1933 vDilepton = vLepton[0]+vLepton[1];
1934 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1935 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1937 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1938 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
1940 if(fChannel == -1) if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) ((TH1D*)(fListPsi2sLoose->At(i)))->Fill(vCandidate.M());
1941 if(fChannel == 1) if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) ((TH1D*)(fListPsi2sLoose->At(i)))->Fill(vCandidate.M());
1946 for(Int_t i=0; i<4; i++){
1947 //cout<<"Tight systematics psi2s, cut"<<i<<endl;
1948 for(Int_t j=0; j<4; j++){
1949 if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
1950 else fJPsiSels[j] = fJPsiSelsMid[j];
1954 nGoodTracks = 0; nSpdHits = 0;
1955 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1956 AliAODTrack *trk = aod->GetTrack(itr);
1957 if( !trk ) continue;
1958 if(!(trk->TestFilterBit(1<<0))) continue;
1960 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1961 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1962 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1963 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1964 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1965 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1967 Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1969 if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1970 if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1971 if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1972 if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1973 if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1975 TrackIndex[nGoodTracks] = itr;
1976 TrackPt[nGoodTracks] = trk->Pt();
1979 if(nGoodTracks > 4) break;
1982 Int_t mass[3]={-1,-1,-1};
1984 nLepton=0; nPion=0; nHighPt=0;
1986 if(nGoodTracks == 4){
1987 if(nSpdHits<2) continue;
1988 MeanPt = GetMedian(TrackPt);
1989 for(Int_t k=0; k<4; k++){
1990 AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);
1991 if(trk->Pt() > MeanPt){
1992 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1993 qLepton[nLepton] = trk->Charge();
1994 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1995 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1998 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1999 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
2005 qPion[nPion] = trk->Charge();
2006 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
2010 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0) && mass[0] != -1 && mass[1] != -1){
2011 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
2012 vDilepton = vLepton[0]+vLepton[1];
2013 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
2014 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
2016 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
2017 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
2019 if(fChannel == -1) if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) ((TH1D*)(fListPsi2sTight->At(i)))->Fill(vCandidate.M());
2020 if(fChannel == 1) if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) ((TH1D*)(fListPsi2sTight->At(i)))->Fill(vCandidate.M());