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