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