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