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