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