Adding Test bit 0 in AOD
[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
54 // my headers
55 #include "AliAnalysisTaskUpcPsi2s.h"
56
57 ClassImp(AliAnalysisTaskUpcPsi2s);
58
59 using std::cout;
60 using std::endl;
61
62 //trees for UPC analysis,
63 // michal.broz@cern.ch
64
65 //_____________________________________________________________________________
66 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s() 
67   : AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fJPsiTree(0),fPsi2sTree(0),
68     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
69     fTOFtrig1(0), fTOFtrig2(0),
70     fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
71     fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
72     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
73     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
74     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
75     fListTrig(0),fHistUpcTriggersPerRun(0),fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0),
76     fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
77     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
78     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
79
80 {
81
82 //Dummy constructor
83
84 }//AliAnalysisTaskUpcPsi2s
85
86
87 //_____________________________________________________________________________
88 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name) 
89   : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fJPsiTree(0),fPsi2sTree(0),
90     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
91     fTOFtrig1(0), fTOFtrig2(0),
92     fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
93     fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
94     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
95     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
96     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
97     fListTrig(0),fHistUpcTriggersPerRun(0),fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0),
98     fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
99     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
100     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
101
102 {
103
104   // Constructor
105   if( strstr(name,"ESD") ) fType = 0;
106   if( strstr(name,"AOD") ) fType = 1;
107   
108   Init();
109
110   DefineOutput(1, TTree::Class());
111   DefineOutput(2, TTree::Class());
112   DefineOutput(3, TList::Class());
113   DefineOutput(4, TList::Class());
114
115 }//AliAnalysisTaskUpcPsi2s
116
117 //_____________________________________________________________________________
118 void AliAnalysisTaskUpcPsi2s::Init()
119 {
120   
121   for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
122   for(Int_t i=0; i<4; i++) fTOFphi[i] = -666;
123   for(Int_t i=0; i<3; i++){
124         fVtxPos[i] = -666; 
125         fVtxErr[i] = -666;
126         fKfVtxPos[i] = -666;
127         }
128
129 }//Init
130
131 //_____________________________________________________________________________
132 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s() 
133 {
134   // Destructor
135   if(fJPsiTree){
136      delete fJPsiTree;
137      fJPsiTree = 0x0;
138   }
139   if(fPsi2sTree){
140      delete fPsi2sTree;
141      fPsi2sTree = 0x0;
142   }
143   if(fListTrig){
144      delete fListTrig;
145      fListTrig = 0x0;
146   }
147   if(fListHist){
148      delete fListHist;
149      fListHist = 0x0;
150   }
151
152 }//~AliAnalysisTaskUpcPsi2s
153
154
155 //_____________________________________________________________________________
156 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
157 {
158   //input file
159   fDataFilnam = new TObjString();
160   fDataFilnam->SetString("");
161
162     //tracks
163   fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
164   fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
165   fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
166   fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
167
168   //output tree with JPsi candidate events
169   fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
170   fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
171   fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
172   fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
173   
174   fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
175   fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
176   fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
177   fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
178   fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
179   fJPsiTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
180   fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
181   
182   fJPsiTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
183   fJPsiTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
184   fJPsiTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
185   
186   fJPsiTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
187   fJPsiTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
188   fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
189   fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
190   
191   fJPsiTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
192   
193   fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
194   fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
195   fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
196   fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");  
197   fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
198   fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
199   fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");                        
200   if( fType == 0 ) {
201     fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
202   }
203   if( fType == 1 ) {
204     fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
205   }
206  
207  //output tree with Psi2s candidate events
208   fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
209   fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
210   fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
211   fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
212   
213   fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
214   fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
215   fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
216   fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
217   fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
218   fPsi2sTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
219   fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
220   
221   fPsi2sTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
222   fPsi2sTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
223   fPsi2sTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
224   
225   fPsi2sTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
226   fPsi2sTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
227   fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
228   fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
229   
230   fPsi2sTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
231   
232   fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
233   fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
234   fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
235   fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");  
236   fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
237   fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
238   fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");                       
239   if( fType == 0 ) {
240     fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
241   }
242   if( fType == 1 ) {
243     fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
244   }
245   
246   fListTrig = new TList();
247   fListTrig ->SetOwner();
248   
249   fHistUpcTriggersPerRun = new TH1D("fHistUpcTriggersPerRun", "fHistUpcTriggersPerRun", 3000, 167000.5, 170000.5);
250   fListTrig->Add(fHistUpcTriggersPerRun);
251   
252   fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 3000, 167000.5, 170000.5);
253   fListTrig->Add(fHistZedTriggersPerRun);
254
255   fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 3000, 167000.5, 170000.5);
256   fListTrig->Add(fHistCvlnTriggersPerRun);
257   
258   fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 3000, 167000.5, 170000.5);
259   fListTrig->Add(fHistMBTriggersPerRun);
260   
261   fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 3000, 167000.5, 170000.5);
262   fListTrig->Add(fHistCentralTriggersPerRun);
263   
264   fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 3000, 167000.5, 170000.5);
265   fListTrig->Add(fHistSemiCentralTriggersPerRun);
266   
267   fListHist = new TList();
268   fListHist ->SetOwner();
269   
270   TString CutNameJPsi[12] = {"Analyzed","Triggered","Vertex cut","V0 decision","Two good tracks",
271                                 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
272   fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",12,0.5,12.5);
273   for (Int_t i = 0; i<12; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
274   fListHist->Add(fHistNeventsJPsi);
275   
276   fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
277   fListHist->Add(fHistTPCsignalJPsi);
278   
279   fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
280   fListHist->Add(fHistDiLeptonPtJPsi);
281
282   fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
283   fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
284   fListHist->Add(fHistDiElectronMass);
285   
286   fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
287   fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
288   fListHist->Add(fHistDiMuonMass);
289
290   TString CutNamePsi2s[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Four good tracks",
291                                 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
292
293   fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",13,0.5,13.5);
294   for (Int_t i = 0; i<13; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
295   fListHist->Add(fHistNeventsPsi2s);
296
297   fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
298   fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
299   fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
300   fListHist->Add(fHistPsi2sMassVsPt);
301   
302   fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
303   fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
304   fListHist->Add(fHistPsi2sMassCoherent);
305   
306   PostData(1, fJPsiTree);
307   PostData(2, fPsi2sTree);
308   PostData(3, fListTrig);
309   PostData(4, fListHist);
310
311 }//UserCreateOutputObjects
312
313 //_____________________________________________________________________________
314 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *) 
315 {
316
317   //cout<<"#################### Next event ##################"<<endl;
318
319   if( fType == 0 ){
320         RunESDtrig(); 
321         if(fRunHist) RunESDhist();
322         if(fRunTree) RunESDtree();
323         }
324
325   if( fType == 1 ){
326         RunAODtrig(); 
327         if(fRunHist) RunAODhist();
328         if(fRunTree) RunAODtree();
329         }
330
331 }//UserExec
332 //_____________________________________________________________________________
333 void AliAnalysisTaskUpcPsi2s::RunAODtrig()
334 {
335
336   //input event
337   AliAODEvent *aod = (AliAODEvent*) InputEvent();
338   if(!aod) return;
339
340   fRunNum = aod ->GetRunNumber();
341   //Trigger
342   TString trigger = aod->GetFiredTriggerClasses();
343   
344   if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
345   
346   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
347   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
348   
349   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
350   if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
351   
352   //MB, Central and SemiCentral triggers
353   AliCentrality *centrality = aod->GetCentrality();
354   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
355   
356   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
357   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
358   
359   if( ((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) || 
360       ((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) || 
361       ((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) ){
362       
363       if(percentile<80 && percentile>0) fHistMBTriggersPerRun->Fill(fRunNum);
364       if(percentile<6 && percentile>0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
365       if(percentile<50 && percentile>15 && (trigger.Contains("CVLN_B2-B"))) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
366       if(percentile<50 && percentile>15 && (trigger.Contains("CVLN_R1-B"))) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
367       
368       }
369     
370 PostData(3, fListTrig);
371
372 }
373 //_____________________________________________________________________________
374 void AliAnalysisTaskUpcPsi2s::RunAODhist()
375 {
376
377   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
378   
379   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
380   Double_t muonMass = partMuon->Mass();
381   
382   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
383   Double_t electronMass = partElectron->Mass();
384   
385   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
386   Double_t pionMass = partPion->Mass();
387
388   //input event
389   AliAODEvent *aod = (AliAODEvent*) InputEvent();
390   if(!aod) return;
391
392   fHistNeventsJPsi->Fill(1);
393   fHistNeventsPsi2s->Fill(1);
394
395   //Trigger
396   TString trigger = aod->GetFiredTriggerClasses();
397   
398   if( !trigger.Contains("CCUP") ) return;
399   
400   fHistNeventsJPsi->Fill(2);
401   fHistNeventsPsi2s->Fill(2);
402
403   //primary vertex
404   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
405   fVtxContrib = fAODVertex->GetNContributors();
406   if(fVtxContrib < 2) return;
407   
408   fHistNeventsJPsi->Fill(3);
409   fHistNeventsPsi2s->Fill(3);
410
411   //VZERO, ZDC
412   AliAODVZERO *fV0data = aod ->GetVZEROData();
413   AliAODZDC *fZDCdata = aod->GetZDCData();
414   
415   fV0Adecision = fV0data->GetV0ADecision();
416   fV0Cdecision = fV0data->GetV0CDecision();
417   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
418   
419   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
420   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
421
422   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
423   
424   fHistNeventsJPsi->Fill(4);
425   fHistNeventsPsi2s->Fill(4);
426
427   //Two tracks loop
428   Int_t nGoodTracks = 0;
429   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
430   
431   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
432   Short_t qLepton[4], qPion[4];
433   UInt_t nLepton=0, nPion=0, nHighPt=0;
434   Double_t fRecTPCsignal[5];
435   Int_t mass[3]={-1,-1,-1};
436   
437    
438   //Four track loop
439   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
440     AliAODTrack *trk = aod->GetTrack(itr);
441     if( !trk ) continue;
442
443       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
444       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
445       if(trk->GetTPCNcls() < 50)continue;
446       if(trk->Chi2perNDF() > 4)continue;
447       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
448       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
449       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
450       delete trk_clone;
451
452       if(TMath::Abs(dca[1]) > 2) continue;
453      
454       TrackIndex[nGoodTracks] = itr;
455       nGoodTracks++;
456                                   
457       if(nGoodTracks > 4) break;  
458   }//Track loop
459   
460   nLepton=0; nPion=0; nHighPt=0;
461   mass[0]= -1; mass[1]= -1, mass[2]= -1;
462   
463   if(nGoodTracks == 4){
464           fHistNeventsPsi2s->Fill(5);
465           for(Int_t i=0; i<4; i++){
466                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
467                 
468                 if(trk->Pt() > 1){   
469                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();      
470                         qLepton[nLepton] = trk->Charge();
471                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
472                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
473                                         mass[nLepton] = 0;
474                                         }
475                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
476                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
477                                         mass[nLepton] = 1;
478                                         }
479                         nLepton++;
480                         }
481                 else{
482                         qPion[nPion] = trk->Charge();
483                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
484                         nPion++;
485                         }
486                                       
487                 if(nLepton > 2 || nPion > 2) break;
488                 }
489         if((nLepton == 2) && (nPion == 2)){
490                 fHistNeventsPsi2s->Fill(6);
491                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
492                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
493                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
494                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
495                         fHistNeventsPsi2s->Fill(10);
496                         if(mass[0] == mass[1]) {
497                                 fHistNeventsPsi2s->Fill(11); 
498                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
499                                 vDilepton = vLepton[0]+vLepton[1];
500                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
501                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
502                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);   
503                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
504                                 }
505                         }
506                 }
507   }
508   
509   nGoodTracks = 0;
510   //Two track loop
511   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
512     AliAODTrack *trk = aod->GetTrack(itr);
513     if( !trk ) continue;
514
515       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
516       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
517       if(trk->GetTPCNcls() < 70)continue;
518       if(trk->Chi2perNDF() > 4)continue;
519       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
520       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
521       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
522       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
523       delete trk_clone;
524       if(TMath::Abs(dca[1]) > 2) continue;
525       if(TMath::Abs(dca[0]) > 0.2) continue;
526      
527       TrackIndex[nGoodTracks] = itr;
528       nGoodTracks++;
529                                   
530       if(nGoodTracks > 2) break;  
531   }//Track loop
532   
533    nLepton=0; nPion=0; nHighPt=0;
534   mass[0]= -1; mass[1]= -1, mass[2]= -1;
535
536   if(nGoodTracks == 2){
537           fHistNeventsJPsi->Fill(5);
538           for(Int_t i=0; i<2; i++){
539                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);                
540                 if(trk->Pt() > 1) nHighPt++;     
541                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();     
542                 qLepton[nLepton] = trk->Charge();
543                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
544                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
545                                 mass[nLepton] = 0;
546                                 }
547                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
548                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
549                                 mass[nLepton] = 1;
550                                 }
551                 nLepton++;              
552                 }               
553         if(nLepton == 2){
554                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
555                 if(qLepton[0]*qLepton[1] < 0){
556                         fHistNeventsJPsi->Fill(7);
557                         if(nHighPt > 0){
558                                 fHistNeventsJPsi->Fill(8);
559                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
560                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
561                                 if(mass[0] == mass[1] && mass[0] != -1) {
562                                         fHistNeventsJPsi->Fill(10);
563                                         vCandidate = vLepton[0]+vLepton[1];
564                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
565                                         if(mass[0] == 0) {
566                                                 fHistDiMuonMass->Fill(vCandidate.M());
567                                                 fHistNeventsJPsi->Fill(11);
568                                                 }
569                                         if(mass[0] == 1) {
570                                                 fHistDiElectronMass->Fill(vCandidate.M());
571                                                 fHistNeventsJPsi->Fill(12);
572                                                 }
573                                         }
574                                 }
575                         }
576                 }
577   }
578   
579  
580   PostData(4, fListHist);
581
582 }
583
584 //_____________________________________________________________________________
585 void AliAnalysisTaskUpcPsi2s::RunAODtree()
586 {
587   //input event
588   AliAODEvent *aod = (AliAODEvent*) InputEvent();
589   if(!aod) return;
590
591   //input data
592   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
593   fDataFilnam->Clear();
594   fDataFilnam->SetString(filnam);
595   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
596   fRunNum = aod ->GetRunNumber();
597
598   //Trigger
599   TString trigger = aod->GetFiredTriggerClasses();
600   
601   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
602   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
603   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
604   
605   Bool_t isTriggered = kFALSE;
606   for(Int_t i=0; i<ntrg; i++) {
607     if( fTrigger[i] ) isTriggered = kTRUE;
608   }
609   if( !isTriggered ) return;
610
611   //trigger inputs
612   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
613   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();  
614
615   //Event identification
616   fPerNum = aod ->GetPeriodNumber();
617   fOrbNum = aod ->GetOrbitNumber();
618   fBCrossNum = aod ->GetBunchCrossNumber();
619
620   //primary vertex
621   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
622   fVtxContrib = fAODVertex->GetNContributors();
623   fVtxPos[0] = fAODVertex->GetX();
624   fVtxPos[1] = fAODVertex->GetY();
625   fVtxPos[2] = fAODVertex->GetZ();
626   Double_t CovMatx[6];
627   fAODVertex->GetCovarianceMatrix(CovMatx); 
628   fVtxErr[0] = CovMatx[0];
629   fVtxErr[1] = CovMatx[1];
630   fVtxErr[2] = CovMatx[2];
631   fVtxChi2 = fAODVertex->GetChi2();
632   fVtxNDF = fAODVertex->GetNDF();
633
634   //Tracklets
635   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
636
637   //VZERO, ZDC
638   AliAODVZERO *fV0data = aod ->GetVZEROData();
639   AliAODZDC *fZDCdata = aod->GetZDCData();
640   
641   fV0Adecision = fV0data->GetV0ADecision();
642   fV0Cdecision = fV0data->GetV0CDecision();
643   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
644   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
645   
646   fNLooseTracks = 0;
647   
648   //Track loop - loose cuts
649   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
650     AliAODTrack *trk = aod->GetTrack(itr);
651     if( !trk ) continue;
652
653       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
654       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
655       if(trk->GetTPCNcls() < 20)continue;
656       fNLooseTracks++; 
657   }//Track loop -loose cuts
658   
659   Int_t nGoodTracks=0;
660   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
661   
662   //Two track loop
663   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
664     AliAODTrack *trk = aod->GetTrack(itr);
665     if( !trk ) continue;
666     if(!(trk->TestFilterBit(1<<0)))continue;
667
668       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
669       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
670       if(trk->GetTPCNcls() < 70)continue;
671       if(trk->Chi2perNDF() > 4)continue;
672       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
673       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
674       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
675       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
676       delete trk_clone;
677       if(TMath::Abs(dca[1]) > 2) continue;
678       if(TMath::Abs(dca[0]) > 0.2) continue;
679      
680       TrackIndex[nGoodTracks] = itr;
681       nGoodTracks++;
682                                   
683       if(nGoodTracks > 2) break;  
684   }//Track loop
685
686   if(nGoodTracks == 2){
687   
688           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
689           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
690           Double_t muonMass = partMuon->Mass();  
691           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
692           Double_t electronMass = partElectron->Mass();  
693           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
694           Double_t pionMass = partPion->Mass();
695   
696           Double_t KFcov[21];
697           Double_t KFpar[6];
698           Double_t KFmass = pionMass;
699           Double_t fRecTPCsignal;
700           AliKFParticle *KFpart[2];
701           AliKFVertex *KFvtx = new AliKFVertex();
702           KFvtx->SetField(aod->GetMagneticField()); 
703   
704           for(Int_t i=0; i<2; i++){
705                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
706                 
707                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
708                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
709                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
710                 delete trk_clone;
711                                 
712                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk); 
713                 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
714                 
715                 trk->GetPosition(KFpar);
716                 trk->PxPyPz(KFpar+3);
717                 trk->GetCovMatrix(KFcov);
718                 
719                 if(trk->Pt() > 1){   
720                         fRecTPCsignal = trk->GetTPCsignal();      
721                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
722                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
723                         }
724                 else KFmass = pionMass;
725                 
726                 KFpart[i] = new AliKFParticle();
727                 KFpart[i]->SetField(aod->GetMagneticField());
728                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
729                 KFvtx->AddDaughter(*KFpart[i]); 
730                 
731                 
732                 Double_t pos[3]={0,0,0};
733                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
734                 parTrk->CopyFromVTrack((AliVTrack*) trk);
735                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
736                 else {
737                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
738                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
739                      }
740                 delete parTrk;          
741                 }
742   fKfVtxPos[0]= KFvtx->GetX();
743   fKfVtxPos[1]= KFvtx->GetY();
744   fKfVtxPos[2]= KFvtx->GetZ();
745   for(UInt_t i=0; i<2; i++)delete KFpart[i];
746   delete KFvtx; 
747
748   fJPsiTree ->Fill();
749   }
750   
751    nGoodTracks = 0;
752    //Four track loop
753   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
754     AliAODTrack *trk = aod->GetTrack(itr);
755     if( !trk ) continue;
756     if(!(trk->TestFilterBit(1<<0)))continue;
757
758       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
759       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
760       if(trk->GetTPCNcls() < 50)continue;
761       if(trk->Chi2perNDF() > 4)continue;
762       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
763       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
764       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
765       delete trk_clone;
766       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
767       if(TMath::Abs(dca[1]) > 2) continue;
768
769       TrackIndex[nGoodTracks] = itr;
770       nGoodTracks++;
771                                   
772       if(nGoodTracks > 4) break;  
773   }//Track loop
774       
775     
776   if(nGoodTracks == 4){
777   
778           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
779           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
780           Double_t muonMass = partMuon->Mass();  
781           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
782           Double_t electronMass = partElectron->Mass();  
783           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
784           Double_t pionMass = partPion->Mass();
785   
786           Double_t KFcov[21];
787           Double_t KFpar[6];
788           Double_t KFmass = pionMass;
789           Double_t fRecTPCsignal;
790           AliKFParticle *KFpart[4];
791           AliKFVertex *KFvtx = new AliKFVertex();
792           KFvtx->SetField(aod->GetMagneticField()); 
793                   
794           for(Int_t i=0; i<4; i++){
795                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
796                 
797                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
798                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
799                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
800                 delete trk_clone;
801                 
802                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
803                 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
804                                 
805                 trk->GetPosition(KFpar);
806                 trk->PxPyPz(KFpar+3);
807                 trk->GetCovMatrix(KFcov);
808                 
809                 if(trk->Pt() > 1){   
810                         fRecTPCsignal = trk->GetTPCsignal();      
811                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
812                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
813                         }
814                 else KFmass = pionMass;
815                 
816                 KFpart[i] = new AliKFParticle();
817                 KFpart[i]->SetField(aod->GetMagneticField());
818                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
819                 KFvtx->AddDaughter(*KFpart[i]); 
820                                 
821                 Double_t pos[3]={0,0,0};
822                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
823                 parTrk->CopyFromVTrack((AliVTrack*) trk);
824                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
825                 else {
826                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
827                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
828                      }
829                 delete parTrk;          
830                 }
831   fKfVtxPos[0]= KFvtx->GetX();
832   fKfVtxPos[1]= KFvtx->GetY();
833   fKfVtxPos[2]= KFvtx->GetZ();
834   for(UInt_t i=0; i<4; i++)delete KFpart[i];
835   delete KFvtx; 
836   fPsi2sTree ->Fill();
837   }
838   
839   PostData(1, fJPsiTree);
840   PostData(2, fPsi2sTree);
841
842 }//RunAOD
843
844 //_____________________________________________________________________________
845 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
846 {
847
848   //input event
849   AliESDEvent *esd = (AliESDEvent*) InputEvent();
850   if(!esd) return;
851
852   fRunNum = esd ->GetRunNumber();
853   //Trigger
854   TString trigger = esd->GetFiredTriggerClasses();
855   
856   if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
857   
858   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
859   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
860   
861   if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
862   
863    //MB, Central and SemiCentral triggers
864   AliCentrality *centrality = esd->GetCentrality();
865   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
866   
867   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
868   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
869   
870   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<80 && percentile>0) fHistMBTriggersPerRun->Fill(fRunNum);
871   
872   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<6 && percentile>0) fHistCentralTriggersPerRun->Fill(fRunNum);
873
874   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<50 && percentile>15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
875
876   
877 PostData(3, fListTrig);
878
879 }
880 //_____________________________________________________________________________
881 void AliAnalysisTaskUpcPsi2s::RunESDhist()
882 {
883
884   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
885   
886   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
887   Double_t muonMass = partMuon->Mass();
888   
889   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
890   Double_t electronMass = partElectron->Mass();
891   
892   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
893   Double_t pionMass = partPion->Mass();
894
895   //input event
896   AliESDEvent *esd = (AliESDEvent*) InputEvent();
897   if(!esd) return;
898
899   fHistNeventsJPsi->Fill(1);
900   fHistNeventsPsi2s->Fill(1);
901
902   //Trigger
903   TString trigger = esd->GetFiredTriggerClasses();
904   
905   if( !trigger.Contains("CCUP") ) return;
906   
907   fHistNeventsJPsi->Fill(2);
908   fHistNeventsPsi2s->Fill(2);
909
910   //primary vertex
911   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
912   fVtxContrib = fESDVertex->GetNContributors();
913   if(fVtxContrib < 2) return;
914   
915   fHistNeventsJPsi->Fill(3);
916   fHistNeventsPsi2s->Fill(3);
917
918   //VZERO, ZDC
919   AliESDVZERO *fV0data = esd->GetVZEROData();
920   AliESDZDC *fZDCdata = esd->GetESDZDC();
921   
922   fV0Adecision = fV0data->GetV0ADecision();
923   fV0Cdecision = fV0data->GetV0CDecision();
924   if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
925   
926   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
927   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
928   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
929   
930   fHistNeventsJPsi->Fill(4);
931   fHistNeventsPsi2s->Fill(4);
932
933    Int_t nGoodTracks=0;
934   //Two tracks loop
935   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
936   
937   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
938   Short_t qLepton[4], qPion[4];
939   UInt_t nLepton=0, nPion=0, nHighPt=0;
940   Double_t fRecTPCsignal[5];
941   Int_t mass[3]={-1,-1,-1};
942
943  //Two Track loop
944   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
945     AliESDtrack *trk = esd->GetTrack(itr);
946     if( !trk ) continue;
947
948       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
949       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
950       if(trk->GetTPCNcls() < 70)continue;
951       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
952       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
953       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
954       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
955       trk->GetImpactParameters(dca[0],dca[1]);
956       if(TMath::Abs(dca[1]) > 2) continue;
957       if(TMath::Abs(dca[1]) > 0.2) continue;
958       
959       TrackIndex[nGoodTracks] = itr;
960       nGoodTracks++;
961       if(nGoodTracks > 2) break;   
962   }//Track loop
963
964   
965   if(nGoodTracks == 2){
966           fHistNeventsJPsi->Fill(5);
967           for(Int_t i=0; i<2; i++){
968                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);                
969                 if(trk->Pt() > 1) nHighPt++;     
970                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();     
971                 qLepton[nLepton] = trk->Charge();
972                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
973                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
974                                 mass[nLepton] = 0;
975                                 }
976                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
977                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
978                                 mass[nLepton] = 1;
979                                 }
980                 nLepton++;              
981                 }               
982         if(nLepton == 2){
983                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
984                 if(qLepton[0]*qLepton[1] < 0){
985                         fHistNeventsJPsi->Fill(7);
986                         if(nHighPt > 0){
987                                 fHistNeventsJPsi->Fill(8);
988                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
989                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
990                                 if(mass[0] == mass[1] && mass[0] != -1) {
991                                         fHistNeventsJPsi->Fill(10);
992                                         vCandidate = vLepton[0]+vLepton[1];
993                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
994                                         if(mass[0] == 0) {
995                                                 fHistDiMuonMass->Fill(vCandidate.M());
996                                                 fHistNeventsJPsi->Fill(11);
997                                                 }
998                                         if(mass[0] == 1) {
999                                                 fHistDiElectronMass->Fill(vCandidate.M());
1000                                                 fHistNeventsJPsi->Fill(12);
1001                                                 }
1002                                         }
1003                                 }
1004                         }
1005                 }
1006   }
1007   nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1008   mass[0]= -1; mass[1]= -1, mass[2]= -1;
1009   
1010     //Four Track loop
1011   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1012     AliESDtrack *trk = esd->GetTrack(itr);
1013     if( !trk ) continue;
1014
1015       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1016       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1017       if(trk->GetTPCNcls() < 50)continue;
1018       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1019       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1020       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1021       trk->GetImpactParameters(dca[0],dca[1]);
1022       if(TMath::Abs(dca[1]) > 2) continue;
1023       
1024       TrackIndex[nGoodTracks] = itr;
1025       nGoodTracks++;
1026       if(nGoodTracks > 4) break;   
1027   }//Track loop
1028   
1029   if(nGoodTracks == 4){
1030           fHistNeventsPsi2s->Fill(5);
1031           for(Int_t i=0; i<4; i++){
1032                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1033                 
1034                 if(trk->Pt() > 1){   
1035                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();      
1036                         qLepton[nLepton] = trk->Charge();
1037                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1038                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1039                                         mass[nLepton] = 0;
1040                                         }
1041                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1042                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1043                                         mass[nLepton] = 1;
1044                                         }
1045                         nLepton++;
1046                         }
1047                 else{
1048                         qPion[nPion] = trk->Charge();
1049                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1050                         nPion++;
1051                         }
1052                                       
1053                 if(nLepton > 2 || nPion > 2) break;
1054                 }
1055         if((nLepton == 2) && (nPion == 2)){
1056                 fHistNeventsPsi2s->Fill(6);
1057                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1058                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1059                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1060                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1061                         fHistNeventsPsi2s->Fill(10);
1062                         if(mass[0] == mass[1]) {
1063                                 fHistNeventsPsi2s->Fill(11); 
1064                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1065                                 vDilepton = vLepton[0]+vLepton[1];
1066                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1067                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1068                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);   
1069                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1070                                 }
1071                         }
1072                 }
1073   }
1074   
1075   PostData(4, fListHist);
1076
1077 }
1078
1079 //_____________________________________________________________________________
1080 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1081 {
1082
1083   //input event
1084   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1085   if(!esd) return;
1086
1087   //input data
1088   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1089   fDataFilnam->Clear();
1090   fDataFilnam->SetString(filnam);
1091   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1092   fRunNum = esd->GetRunNumber();
1093
1094    //Trigger
1095   TString trigger = esd->GetFiredTriggerClasses();
1096   
1097   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1098   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
1099   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1100   
1101   Bool_t isTriggered = kFALSE;
1102   for(Int_t i=0; i<ntrg; i++) {
1103     if( fTrigger[i] ) isTriggered = kTRUE;
1104   }
1105   if( !isTriggered ) return;
1106   
1107   //trigger inputs
1108   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1109   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1110   
1111   //TOF trigger info (0OMU)
1112   fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1113   fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1114
1115   //Event identification
1116   fPerNum = esd->GetPeriodNumber();
1117   fOrbNum = esd->GetOrbitNumber();
1118   fBCrossNum = esd->GetBunchCrossNumber();
1119
1120   //primary vertex
1121   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1122   fVtxContrib = fESDVertex->GetNContributors();
1123   fVtxPos[0] = fESDVertex->GetX();
1124   fVtxPos[1] = fESDVertex->GetY();
1125   fVtxPos[2] = fESDVertex->GetZ();
1126   Double_t CovMatx[6];
1127   fESDVertex->GetCovarianceMatrix(CovMatx); 
1128   fVtxErr[0] = CovMatx[0];
1129   fVtxErr[1] = CovMatx[1];
1130   fVtxErr[2] = CovMatx[2];
1131   fVtxChi2 = fESDVertex->GetChi2();
1132   fVtxNDF = fESDVertex->GetNDF();
1133
1134   //Tracklets
1135   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1136
1137   //VZERO, ZDC
1138   AliESDVZERO *fV0data = esd->GetVZEROData();
1139   AliESDZDC *fZDCdata = esd->GetESDZDC();
1140   
1141   fV0Adecision = fV0data->GetV0ADecision();
1142   fV0Cdecision = fV0data->GetV0CDecision();
1143   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1144   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1145
1146   fNLooseTracks = 0;
1147   
1148   //Track loop - loose cuts
1149   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1150     AliESDtrack *trk = esd->GetTrack(itr);
1151     if( !trk ) continue;
1152
1153       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1154       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1155       if(trk->GetTPCNcls() < 20)continue;
1156       fNLooseTracks++; 
1157   }//Track loop -loose cuts
1158   
1159   Int_t nGoodTracks=0;
1160   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1161   
1162   //Two Track loop
1163   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1164     AliESDtrack *trk = esd->GetTrack(itr);
1165     if( !trk ) continue;
1166
1167       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1168       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1169       if(trk->GetTPCNcls() < 70)continue;
1170       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1171       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1172       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1173       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1174       trk->GetImpactParameters(dca[0],dca[1]);
1175       if(TMath::Abs(dca[1]) > 2) continue;
1176       if(TMath::Abs(dca[1]) > 0.2) continue;
1177       
1178       TrackIndex[nGoodTracks] = itr;
1179       nGoodTracks++;
1180       if(nGoodTracks > 2) break;   
1181   }//Track loop
1182
1183   if(nGoodTracks == 2){
1184           for(Int_t i=0; i<2; i++){
1185                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1186                 
1187                 AliExternalTrackParam cParam;
1188                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1189                                 
1190                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk); 
1191                 
1192                 Double_t pos[3]={0,0,0};
1193                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1194                 else {
1195                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1196                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1197                      }          
1198                 }
1199   fJPsiTree ->Fill();
1200   }
1201   
1202   nGoodTracks = 0;
1203   //Four track loop
1204   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1205     AliESDtrack *trk = esd->GetTrack(itr);
1206     if( !trk ) continue;
1207
1208       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1209       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1210       if(trk->GetTPCNcls() < 50)continue;
1211       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1212       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1213       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1214       trk->GetImpactParameters(dca[0],dca[1]);
1215       if(TMath::Abs(dca[1]) > 2) continue;
1216       
1217       TrackIndex[nGoodTracks] = itr;
1218       nGoodTracks++;
1219       if(nGoodTracks > 4) break;   
1220   }//Track loop
1221   
1222   if(nGoodTracks == 4){
1223           for(Int_t i=0; i<4; i++){
1224                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1225                 
1226                 AliExternalTrackParam cParam;
1227                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1228
1229                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1230                  
1231                 Double_t pos[3]={0,0,0};
1232                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1233                 else {
1234                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1235                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1236                      }          
1237                 }
1238   fPsi2sTree ->Fill();
1239   }
1240   
1241   PostData(1, fJPsiTree);
1242   PostData(2, fPsi2sTree);
1243
1244 }//RunESD
1245
1246 //_____________________________________________________________________________
1247 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *) 
1248 {
1249
1250   cout<<"Analysis complete."<<endl;
1251 }//Terminate
1252