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