811bd65d89f5bacc55a4d2a1196700e2f8e8a583
[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) && percentile<80 && percentile>0) fHistMBTriggersPerRun->Fill(fRunNum);
360   
361   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<6 && percentile>0) fHistCentralTriggersPerRun->Fill(fRunNum);
362
363   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<50 && percentile>15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
364
365 PostData(3, fListTrig);
366
367 }
368 //_____________________________________________________________________________
369 void AliAnalysisTaskUpcPsi2s::RunAODhist()
370 {
371
372   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
373   
374   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
375   Double_t muonMass = partMuon->Mass();
376   
377   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
378   Double_t electronMass = partElectron->Mass();
379   
380   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
381   Double_t pionMass = partPion->Mass();
382
383   //input event
384   AliAODEvent *aod = (AliAODEvent*) InputEvent();
385   if(!aod) return;
386
387   fHistNeventsJPsi->Fill(1);
388   fHistNeventsPsi2s->Fill(1);
389
390   //Trigger
391   TString trigger = aod->GetFiredTriggerClasses();
392   
393   if( !trigger.Contains("CCUP") ) return;
394   
395   fHistNeventsJPsi->Fill(2);
396   fHistNeventsPsi2s->Fill(2);
397
398   //primary vertex
399   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
400   fVtxContrib = fAODVertex->GetNContributors();
401   if(fVtxContrib < 2) return;
402   
403   fHistNeventsJPsi->Fill(3);
404   fHistNeventsPsi2s->Fill(3);
405
406   //VZERO, ZDC
407   AliAODVZERO *fV0data = aod ->GetVZEROData();
408   AliAODZDC *fZDCdata = aod->GetZDCData();
409   
410   fV0Adecision = fV0data->GetV0ADecision();
411   fV0Cdecision = fV0data->GetV0CDecision();
412   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
413   
414   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
415   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
416
417   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
418   
419   fHistNeventsJPsi->Fill(4);
420   fHistNeventsPsi2s->Fill(4);
421
422   //Two tracks loop
423   Int_t nGoodTracks = 0;
424   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
425   
426   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
427   Short_t qLepton[4], qPion[4];
428   UInt_t nLepton=0, nPion=0, nHighPt=0;
429   Double_t fRecTPCsignal[5];
430   Int_t mass[3]={-1,-1,-1};
431   
432    
433   //Four track loop
434   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
435     AliAODTrack *trk = aod->GetTrack(itr);
436     if( !trk ) continue;
437
438       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
439       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
440       if(trk->GetTPCNcls() < 50)continue;
441       if(trk->Chi2perNDF() > 4)continue;
442       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
443       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
444       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
445       delete trk_clone;
446
447       if(TMath::Abs(dca[1]) > 2) continue;
448      
449       TrackIndex[nGoodTracks] = itr;
450       nGoodTracks++;
451                                   
452       if(nGoodTracks > 4) break;  
453   }//Track loop
454   
455   nLepton=0; nPion=0; nHighPt=0;
456   mass[0]= -1; mass[1]= -1, mass[2]= -1;
457   
458   if(nGoodTracks == 4){
459           fHistNeventsPsi2s->Fill(5);
460           for(Int_t i=0; i<4; i++){
461                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
462                 
463                 if(trk->Pt() > 1){   
464                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();      
465                         qLepton[nLepton] = trk->Charge();
466                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
467                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
468                                         mass[nLepton] = 0;
469                                         }
470                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
471                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
472                                         mass[nLepton] = 1;
473                                         }
474                         nLepton++;
475                         }
476                 else{
477                         qPion[nPion] = trk->Charge();
478                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
479                         nPion++;
480                         }
481                                       
482                 if(nLepton > 2 || nPion > 2) break;
483                 }
484         if((nLepton == 2) && (nPion == 2)){
485                 fHistNeventsPsi2s->Fill(6);
486                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
487                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
488                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
489                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
490                         fHistNeventsPsi2s->Fill(10);
491                         if(mass[0] == mass[1]) {
492                                 fHistNeventsPsi2s->Fill(11); 
493                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
494                                 vDilepton = vLepton[0]+vLepton[1];
495                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
496                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
497                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);   
498                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
499                                 }
500                         }
501                 }
502   }
503   
504   nGoodTracks = 0;
505   //Two track loop
506   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
507     AliAODTrack *trk = aod->GetTrack(itr);
508     if( !trk ) continue;
509
510       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
511       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
512       if(trk->GetTPCNcls() < 70)continue;
513       if(trk->Chi2perNDF() > 4)continue;
514       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
515       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
516       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
517       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
518       delete trk_clone;
519       if(TMath::Abs(dca[1]) > 2) continue;
520       if(TMath::Abs(dca[0]) > 0.2) continue;
521      
522       TrackIndex[nGoodTracks] = itr;
523       nGoodTracks++;
524                                   
525       if(nGoodTracks > 2) break;  
526   }//Track loop
527   
528    nLepton=0; nPion=0; nHighPt=0;
529   mass[0]= -1; mass[1]= -1, mass[2]= -1;
530
531   if(nGoodTracks == 2){
532           fHistNeventsJPsi->Fill(5);
533           for(Int_t i=0; i<2; i++){
534                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);                
535                 if(trk->Pt() > 1) nHighPt++;     
536                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();     
537                 qLepton[nLepton] = trk->Charge();
538                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
539                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
540                                 mass[nLepton] = 0;
541                                 }
542                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
543                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
544                                 mass[nLepton] = 1;
545                                 }
546                 nLepton++;              
547                 }               
548         if(nLepton == 2){
549                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
550                 if(qLepton[0]*qLepton[1] < 0){
551                         fHistNeventsJPsi->Fill(7);
552                         if(nHighPt > 0){
553                                 fHistNeventsJPsi->Fill(8);
554                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
555                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
556                                 if(mass[0] == mass[1] && mass[0] != -1) {
557                                         fHistNeventsJPsi->Fill(10);
558                                         vCandidate = vLepton[0]+vLepton[1];
559                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
560                                         if(mass[0] == 0) {
561                                                 fHistDiMuonMass->Fill(vCandidate.M());
562                                                 fHistNeventsJPsi->Fill(11);
563                                                 }
564                                         if(mass[0] == 1) {
565                                                 fHistDiElectronMass->Fill(vCandidate.M());
566                                                 fHistNeventsJPsi->Fill(12);
567                                                 }
568                                         }
569                                 }
570                         }
571                 }
572   }
573   
574  
575   PostData(4, fListHist);
576
577 }
578
579 //_____________________________________________________________________________
580 void AliAnalysisTaskUpcPsi2s::RunAODtree()
581 {
582   //input event
583   AliAODEvent *aod = (AliAODEvent*) InputEvent();
584   if(!aod) return;
585
586   //input data
587   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
588   fDataFilnam->Clear();
589   fDataFilnam->SetString(filnam);
590   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
591   fRunNum = aod ->GetRunNumber();
592
593   //Trigger
594   TString trigger = aod->GetFiredTriggerClasses();
595   
596   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
597   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
598   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
599   
600   Bool_t isTriggered = kFALSE;
601   for(Int_t i=0; i<ntrg; i++) {
602     if( fTrigger[i] ) isTriggered = kTRUE;
603   }
604   if( !isTriggered ) return;
605
606   //trigger inputs
607   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
608   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();  
609
610   //Event identification
611   fPerNum = aod ->GetPeriodNumber();
612   fOrbNum = aod ->GetOrbitNumber();
613   fBCrossNum = aod ->GetBunchCrossNumber();
614
615   //primary vertex
616   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
617   fVtxContrib = fAODVertex->GetNContributors();
618   fVtxPos[0] = fAODVertex->GetX();
619   fVtxPos[1] = fAODVertex->GetY();
620   fVtxPos[2] = fAODVertex->GetZ();
621   Double_t CovMatx[6];
622   fAODVertex->GetCovarianceMatrix(CovMatx); 
623   fVtxErr[0] = CovMatx[0];
624   fVtxErr[1] = CovMatx[1];
625   fVtxErr[2] = CovMatx[2];
626   fVtxChi2 = fAODVertex->GetChi2();
627   fVtxNDF = fAODVertex->GetNDF();
628
629   //Tracklets
630   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
631
632   //VZERO, ZDC
633   AliAODVZERO *fV0data = aod ->GetVZEROData();
634   AliAODZDC *fZDCdata = aod->GetZDCData();
635   
636   fV0Adecision = fV0data->GetV0ADecision();
637   fV0Cdecision = fV0data->GetV0CDecision();
638   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
639   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
640   
641   fNLooseTracks = 0;
642   
643   //Track loop - loose cuts
644   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
645     AliAODTrack *trk = aod->GetTrack(itr);
646     if( !trk ) continue;
647
648       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
649       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
650       if(trk->GetTPCNcls() < 20)continue;
651       fNLooseTracks++; 
652   }//Track loop -loose cuts
653   
654   Int_t nGoodTracks=0;
655   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
656   
657   //Two track loop
658   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
659     AliAODTrack *trk = aod->GetTrack(itr);
660     if( !trk ) continue;
661
662       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
663       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
664       if(trk->GetTPCNcls() < 70)continue;
665       if(trk->Chi2perNDF() > 4)continue;
666       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
667       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
668       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
669       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
670       delete trk_clone;
671       if(TMath::Abs(dca[1]) > 2) continue;
672       if(TMath::Abs(dca[0]) > 0.2) continue;
673      
674       TrackIndex[nGoodTracks] = itr;
675       nGoodTracks++;
676                                   
677       if(nGoodTracks > 2) break;  
678   }//Track loop
679
680   if(nGoodTracks == 2){
681   
682           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
683           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
684           Double_t muonMass = partMuon->Mass();  
685           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
686           Double_t electronMass = partElectron->Mass();  
687           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
688           Double_t pionMass = partPion->Mass();
689   
690           Double_t KFcov[21];
691           Double_t KFpar[6];
692           Double_t KFmass = pionMass;
693           Double_t fRecTPCsignal;
694           AliKFParticle *KFpart[2];
695           AliKFVertex *KFvtx = new AliKFVertex();
696           KFvtx->SetField(aod->GetMagneticField()); 
697   
698           for(Int_t i=0; i<2; i++){
699                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
700                 
701                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
702                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
703                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
704                 delete trk_clone;
705                                 
706                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk); 
707                 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
708                 
709                 trk->GetPosition(KFpar);
710                 trk->PxPyPz(KFpar+3);
711                 trk->GetCovMatrix(KFcov);
712                 
713                 if(trk->Pt() > 1){   
714                         fRecTPCsignal = trk->GetTPCsignal();      
715                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
716                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
717                         }
718                 else KFmass = pionMass;
719                 
720                 KFpart[i] = new AliKFParticle();
721                 KFpart[i]->SetField(aod->GetMagneticField());
722                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
723                 KFvtx->AddDaughter(*KFpart[i]); 
724                 
725                 
726                 Double_t pos[3]={0,0,0};
727                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
728                 parTrk->CopyFromVTrack((AliVTrack*) trk);
729                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
730                 else {
731                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
732                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
733                      }
734                 delete parTrk;          
735                 }
736   fKfVtxPos[0]= KFvtx->GetX();
737   fKfVtxPos[1]= KFvtx->GetY();
738   fKfVtxPos[2]= KFvtx->GetZ();
739   for(UInt_t i=0; i<2; i++)delete KFpart[i];
740   delete KFvtx; 
741
742   fJPsiTree ->Fill();
743   }
744   
745    nGoodTracks = 0;
746    //Four track loop
747   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
748     AliAODTrack *trk = aod->GetTrack(itr);
749     if( !trk ) continue;
750
751       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
752       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
753       if(trk->GetTPCNcls() < 50)continue;
754       if(trk->Chi2perNDF() > 4)continue;
755       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
756       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
757       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
758       delete trk_clone;
759       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
760       if(TMath::Abs(dca[1]) > 2) continue;
761
762       TrackIndex[nGoodTracks] = itr;
763       nGoodTracks++;
764                                   
765       if(nGoodTracks > 4) break;  
766   }//Track loop
767       
768     
769   if(nGoodTracks == 4){
770   
771           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
772           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
773           Double_t muonMass = partMuon->Mass();  
774           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
775           Double_t electronMass = partElectron->Mass();  
776           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
777           Double_t pionMass = partPion->Mass();
778   
779           Double_t KFcov[21];
780           Double_t KFpar[6];
781           Double_t KFmass = pionMass;
782           Double_t fRecTPCsignal;
783           AliKFParticle *KFpart[4];
784           AliKFVertex *KFvtx = new AliKFVertex();
785           KFvtx->SetField(aod->GetMagneticField()); 
786                   
787           for(Int_t i=0; i<4; i++){
788                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
789                 
790                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
791                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
792                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
793                 delete trk_clone;
794                 
795                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
796                 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
797                                 
798                 trk->GetPosition(KFpar);
799                 trk->PxPyPz(KFpar+3);
800                 trk->GetCovMatrix(KFcov);
801                 
802                 if(trk->Pt() > 1){   
803                         fRecTPCsignal = trk->GetTPCsignal();      
804                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
805                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
806                         }
807                 else KFmass = pionMass;
808                 
809                 KFpart[i] = new AliKFParticle();
810                 KFpart[i]->SetField(aod->GetMagneticField());
811                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
812                 KFvtx->AddDaughter(*KFpart[i]); 
813                                 
814                 Double_t pos[3]={0,0,0};
815                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
816                 parTrk->CopyFromVTrack((AliVTrack*) trk);
817                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
818                 else {
819                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
820                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
821                      }
822                 delete parTrk;          
823                 }
824   fKfVtxPos[0]= KFvtx->GetX();
825   fKfVtxPos[1]= KFvtx->GetY();
826   fKfVtxPos[2]= KFvtx->GetZ();
827   for(UInt_t i=0; i<4; i++)delete KFpart[i];
828   delete KFvtx; 
829   fPsi2sTree ->Fill();
830   }
831   
832   PostData(1, fJPsiTree);
833   PostData(2, fPsi2sTree);
834
835 }//RunAOD
836
837 //_____________________________________________________________________________
838 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
839 {
840
841   //input event
842   AliESDEvent *esd = (AliESDEvent*) InputEvent();
843   if(!esd) return;
844
845   fRunNum = esd ->GetRunNumber();
846   //Trigger
847   TString trigger = esd->GetFiredTriggerClasses();
848   
849   if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
850   
851   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
852   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
853   
854   if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
855   
856    //MB, Central and SemiCentral triggers
857   AliCentrality *centrality = esd->GetCentrality();
858   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
859   
860   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
861   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
862   
863   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<80 && percentile>0) fHistMBTriggersPerRun->Fill(fRunNum);
864   
865   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<6 && percentile>0) fHistCentralTriggersPerRun->Fill(fRunNum);
866
867   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<50 && percentile>15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
868
869   
870 PostData(3, fListTrig);
871
872 }
873 //_____________________________________________________________________________
874 void AliAnalysisTaskUpcPsi2s::RunESDhist()
875 {
876
877   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
878   
879   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
880   Double_t muonMass = partMuon->Mass();
881   
882   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
883   Double_t electronMass = partElectron->Mass();
884   
885   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
886   Double_t pionMass = partPion->Mass();
887
888   //input event
889   AliESDEvent *esd = (AliESDEvent*) InputEvent();
890   if(!esd) return;
891
892   fHistNeventsJPsi->Fill(1);
893   fHistNeventsPsi2s->Fill(1);
894
895   //Trigger
896   TString trigger = esd->GetFiredTriggerClasses();
897   
898   if( !trigger.Contains("CCUP") ) return;
899   
900   fHistNeventsJPsi->Fill(2);
901   fHistNeventsPsi2s->Fill(2);
902
903   //primary vertex
904   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
905   fVtxContrib = fESDVertex->GetNContributors();
906   if(fVtxContrib < 2) return;
907   
908   fHistNeventsJPsi->Fill(3);
909   fHistNeventsPsi2s->Fill(3);
910
911   //VZERO, ZDC
912   AliESDVZERO *fV0data = esd->GetVZEROData();
913   AliESDZDC *fZDCdata = esd->GetESDZDC();
914   
915   fV0Adecision = fV0data->GetV0ADecision();
916   fV0Cdecision = fV0data->GetV0CDecision();
917   if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
918   
919   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
920   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
921   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
922   
923   fHistNeventsJPsi->Fill(4);
924   fHistNeventsPsi2s->Fill(4);
925
926    Int_t nGoodTracks=0;
927   //Two tracks loop
928   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
929   
930   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
931   Short_t qLepton[4], qPion[4];
932   UInt_t nLepton=0, nPion=0, nHighPt=0;
933   Double_t fRecTPCsignal[5];
934   Int_t mass[3]={-1,-1,-1};
935
936  //Two Track loop
937   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
938     AliESDtrack *trk = esd->GetTrack(itr);
939     if( !trk ) continue;
940
941       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
942       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
943       if(trk->GetTPCNcls() < 70)continue;
944       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
945       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
946       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
947       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
948       trk->GetImpactParameters(dca[0],dca[1]);
949       if(TMath::Abs(dca[1]) > 2) continue;
950       if(TMath::Abs(dca[1]) > 0.2) continue;
951       
952       TrackIndex[nGoodTracks] = itr;
953       nGoodTracks++;
954       if(nGoodTracks > 2) break;   
955   }//Track loop
956
957   
958   if(nGoodTracks == 2){
959           fHistNeventsJPsi->Fill(5);
960           for(Int_t i=0; i<2; i++){
961                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);                
962                 if(trk->Pt() > 1) nHighPt++;     
963                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();     
964                 qLepton[nLepton] = trk->Charge();
965                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
966                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
967                                 mass[nLepton] = 0;
968                                 }
969                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
970                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
971                                 mass[nLepton] = 1;
972                                 }
973                 nLepton++;              
974                 }               
975         if(nLepton == 2){
976                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
977                 if(qLepton[0]*qLepton[1] < 0){
978                         fHistNeventsJPsi->Fill(7);
979                         if(nHighPt > 0){
980                                 fHistNeventsJPsi->Fill(8);
981                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
982                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
983                                 if(mass[0] == mass[1] && mass[0] != -1) {
984                                         fHistNeventsJPsi->Fill(10);
985                                         vCandidate = vLepton[0]+vLepton[1];
986                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
987                                         if(mass[0] == 0) {
988                                                 fHistDiMuonMass->Fill(vCandidate.M());
989                                                 fHistNeventsJPsi->Fill(11);
990                                                 }
991                                         if(mass[0] == 1) {
992                                                 fHistDiElectronMass->Fill(vCandidate.M());
993                                                 fHistNeventsJPsi->Fill(12);
994                                                 }
995                                         }
996                                 }
997                         }
998                 }
999   }
1000   nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1001   mass[0]= -1; mass[1]= -1, mass[2]= -1;
1002   
1003     //Four Track loop
1004   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1005     AliESDtrack *trk = esd->GetTrack(itr);
1006     if( !trk ) continue;
1007
1008       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1009       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1010       if(trk->GetTPCNcls() < 50)continue;
1011       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1012       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1013       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1014       trk->GetImpactParameters(dca[0],dca[1]);
1015       if(TMath::Abs(dca[1]) > 2) continue;
1016       
1017       TrackIndex[nGoodTracks] = itr;
1018       nGoodTracks++;
1019       if(nGoodTracks > 4) break;   
1020   }//Track loop
1021   
1022   if(nGoodTracks == 4){
1023           fHistNeventsPsi2s->Fill(5);
1024           for(Int_t i=0; i<4; i++){
1025                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1026                 
1027                 if(trk->Pt() > 1){   
1028                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();      
1029                         qLepton[nLepton] = trk->Charge();
1030                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1031                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1032                                         mass[nLepton] = 0;
1033                                         }
1034                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1035                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1036                                         mass[nLepton] = 1;
1037                                         }
1038                         nLepton++;
1039                         }
1040                 else{
1041                         qPion[nPion] = trk->Charge();
1042                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1043                         nPion++;
1044                         }
1045                                       
1046                 if(nLepton > 2 || nPion > 2) break;
1047                 }
1048         if((nLepton == 2) && (nPion == 2)){
1049                 fHistNeventsPsi2s->Fill(6);
1050                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1051                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1052                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1053                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1054                         fHistNeventsPsi2s->Fill(10);
1055                         if(mass[0] == mass[1]) {
1056                                 fHistNeventsPsi2s->Fill(11); 
1057                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1058                                 vDilepton = vLepton[0]+vLepton[1];
1059                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1060                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1061                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);   
1062                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1063                                 }
1064                         }
1065                 }
1066   }
1067   
1068   PostData(4, fListHist);
1069
1070 }
1071
1072 //_____________________________________________________________________________
1073 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1074 {
1075
1076   //input event
1077   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1078   if(!esd) return;
1079
1080   //input data
1081   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1082   fDataFilnam->Clear();
1083   fDataFilnam->SetString(filnam);
1084   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1085   fRunNum = esd->GetRunNumber();
1086
1087    //Trigger
1088   TString trigger = esd->GetFiredTriggerClasses();
1089   
1090   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1091   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
1092   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1093   
1094   Bool_t isTriggered = kFALSE;
1095   for(Int_t i=0; i<ntrg; i++) {
1096     if( fTrigger[i] ) isTriggered = kTRUE;
1097   }
1098   if( !isTriggered ) return;
1099   
1100   //trigger inputs
1101   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1102   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1103   
1104   //TOF trigger info (0OMU)
1105   fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1106   fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1107
1108   //Event identification
1109   fPerNum = esd->GetPeriodNumber();
1110   fOrbNum = esd->GetOrbitNumber();
1111   fBCrossNum = esd->GetBunchCrossNumber();
1112
1113   //primary vertex
1114   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1115   fVtxContrib = fESDVertex->GetNContributors();
1116   fVtxPos[0] = fESDVertex->GetX();
1117   fVtxPos[1] = fESDVertex->GetY();
1118   fVtxPos[2] = fESDVertex->GetZ();
1119   Double_t CovMatx[6];
1120   fESDVertex->GetCovarianceMatrix(CovMatx); 
1121   fVtxErr[0] = CovMatx[0];
1122   fVtxErr[1] = CovMatx[1];
1123   fVtxErr[2] = CovMatx[2];
1124   fVtxChi2 = fESDVertex->GetChi2();
1125   fVtxNDF = fESDVertex->GetNDF();
1126
1127   //Tracklets
1128   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1129
1130   //VZERO, ZDC
1131   AliESDVZERO *fV0data = esd->GetVZEROData();
1132   AliESDZDC *fZDCdata = esd->GetESDZDC();
1133   
1134   fV0Adecision = fV0data->GetV0ADecision();
1135   fV0Cdecision = fV0data->GetV0CDecision();
1136   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1137   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1138
1139   fNLooseTracks = 0;
1140   
1141   //Track loop - loose cuts
1142   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1143     AliESDtrack *trk = esd->GetTrack(itr);
1144     if( !trk ) continue;
1145
1146       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1147       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1148       if(trk->GetTPCNcls() < 20)continue;
1149       fNLooseTracks++; 
1150   }//Track loop -loose cuts
1151   
1152   Int_t nGoodTracks=0;
1153   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1154   
1155   //Two Track loop
1156   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1157     AliESDtrack *trk = esd->GetTrack(itr);
1158     if( !trk ) continue;
1159
1160       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1161       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1162       if(trk->GetTPCNcls() < 70)continue;
1163       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1164       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1165       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1166       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1167       trk->GetImpactParameters(dca[0],dca[1]);
1168       if(TMath::Abs(dca[1]) > 2) continue;
1169       if(TMath::Abs(dca[1]) > 0.2) continue;
1170       
1171       TrackIndex[nGoodTracks] = itr;
1172       nGoodTracks++;
1173       if(nGoodTracks > 2) break;   
1174   }//Track loop
1175
1176   if(nGoodTracks == 2){
1177           for(Int_t i=0; i<2; i++){
1178                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1179                 
1180                 AliExternalTrackParam cParam;
1181                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1182                                 
1183                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk); 
1184                 
1185                 Double_t pos[3]={0,0,0};
1186                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1187                 else {
1188                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1189                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1190                      }          
1191                 }
1192   fJPsiTree ->Fill();
1193   }
1194   
1195   nGoodTracks = 0;
1196   //Four track loop
1197   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1198     AliESDtrack *trk = esd->GetTrack(itr);
1199     if( !trk ) continue;
1200
1201       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1202       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1203       if(trk->GetTPCNcls() < 50)continue;
1204       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1205       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1206       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1207       trk->GetImpactParameters(dca[0],dca[1]);
1208       if(TMath::Abs(dca[1]) > 2) continue;
1209       
1210       TrackIndex[nGoodTracks] = itr;
1211       nGoodTracks++;
1212       if(nGoodTracks > 4) break;   
1213   }//Track loop
1214   
1215   if(nGoodTracks == 4){
1216           for(Int_t i=0; i<4; i++){
1217                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1218                 
1219                 AliExternalTrackParam cParam;
1220                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1221
1222                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1223                  
1224                 Double_t pos[3]={0,0,0};
1225                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1226                 else {
1227                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1228                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1229                      }          
1230                 }
1231   fPsi2sTree ->Fill();
1232   }
1233   
1234   PostData(1, fJPsiTree);
1235   PostData(2, fPsi2sTree);
1236
1237 }//RunESD
1238
1239 //_____________________________________________________________________________
1240 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *) 
1241 {
1242
1243   cout<<"Analysis complete."<<endl;
1244 }//Terminate
1245