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