1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
23 #include "TClonesArray.h"
24 #include "TParticle.h"
25 #include "TObjString.h"
27 #include "TDatabasePDG.h"
28 #include "TLorentzVector.h"
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 #include "AliKFVertex.h"
52 #include "AliExternalTrackParam.h"
55 #include "AliAnalysisTaskUpcPsi2s.h"
57 ClassImp(AliAnalysisTaskUpcPsi2s);
62 //trees for UPC analysis,
63 // michal.broz@cern.ch
65 //_____________________________________________________________________________
66 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s()
67 : AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fJPsiTree(0),fPsi2sTree(0),
68 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
69 fTOFtrig1(0), fTOFtrig2(0),
70 fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
71 fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
72 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
73 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
74 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
75 fListTrig(0),fHistUpcTriggersPerRun(0),fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0),
76 fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
77 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
78 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
84 }//AliAnalysisTaskUpcPsi2s
87 //_____________________________________________________________________________
88 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
89 : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fJPsiTree(0),fPsi2sTree(0),
90 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
91 fTOFtrig1(0), fTOFtrig2(0),
92 fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
93 fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
94 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
95 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
96 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
97 fListTrig(0),fHistUpcTriggersPerRun(0),fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0),
98 fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
99 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
100 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
105 if( strstr(name,"ESD") ) fType = 0;
106 if( strstr(name,"AOD") ) fType = 1;
110 DefineOutput(1, TTree::Class());
111 DefineOutput(2, TTree::Class());
112 DefineOutput(3, TList::Class());
113 DefineOutput(4, TList::Class());
115 }//AliAnalysisTaskUpcPsi2s
117 //_____________________________________________________________________________
118 void AliAnalysisTaskUpcPsi2s::Init()
121 for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
122 for(Int_t i=0; i<4; i++) fTOFphi[i] = -666;
123 for(Int_t i=0; i<3; i++){
131 //_____________________________________________________________________________
132 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
152 }//~AliAnalysisTaskUpcPsi2s
155 //_____________________________________________________________________________
156 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
159 fDataFilnam = new TObjString();
160 fDataFilnam->SetString("");
163 fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
164 fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
165 fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
166 fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
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");
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");
179 fJPsiTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
180 fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
182 fJPsiTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
183 fJPsiTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
184 fJPsiTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
186 fJPsiTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
187 fJPsiTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
188 fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
189 fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
191 fJPsiTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
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");
201 fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
204 fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
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");
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");
218 fPsi2sTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
219 fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
221 fPsi2sTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
222 fPsi2sTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
223 fPsi2sTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
225 fPsi2sTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
226 fPsi2sTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
227 fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
228 fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
230 fPsi2sTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
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");
240 fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
243 fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
246 fListTrig = new TList();
247 fListTrig ->SetOwner();
249 fHistUpcTriggersPerRun = new TH1D("fHistUpcTriggersPerRun", "fHistUpcTriggersPerRun", 3000, 167000.5, 170000.5);
250 fListTrig->Add(fHistUpcTriggersPerRun);
252 fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 3000, 167000.5, 170000.5);
253 fListTrig->Add(fHistZedTriggersPerRun);
255 fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 3000, 167000.5, 170000.5);
256 fListTrig->Add(fHistCvlnTriggersPerRun);
258 fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 3000, 167000.5, 170000.5);
259 fListTrig->Add(fHistMBTriggersPerRun);
261 fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 3000, 167000.5, 170000.5);
262 fListTrig->Add(fHistCentralTriggersPerRun);
264 fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 3000, 167000.5, 170000.5);
265 fListTrig->Add(fHistSemiCentralTriggersPerRun);
267 fListHist = new TList();
268 fListHist ->SetOwner();
270 TString CutNameJPsi[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Two good tracks",
271 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
272 fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",13,0.5,13.5);
273 for (Int_t i = 0; i<13; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
274 fListHist->Add(fHistNeventsJPsi);
276 fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
277 fListHist->Add(fHistTPCsignalJPsi);
279 fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
280 fListHist->Add(fHistDiLeptonPtJPsi);
282 fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
283 fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
284 fListHist->Add(fHistDiElectronMass);
286 fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
287 fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
288 fListHist->Add(fHistDiMuonMass);
290 TString CutNamePsi2s[14] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Four good tracks",
291 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
293 fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",14,0.5,14.5);
294 for (Int_t i = 0; i<14; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
295 fListHist->Add(fHistNeventsPsi2s);
297 fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
298 fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
299 fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
300 fListHist->Add(fHistPsi2sMassVsPt);
302 fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
303 fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
304 fListHist->Add(fHistPsi2sMassCoherent);
306 PostData(1, fJPsiTree);
307 PostData(2, fPsi2sTree);
308 PostData(3, fListTrig);
309 PostData(4, fListHist);
311 }//UserCreateOutputObjects
313 //_____________________________________________________________________________
314 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
317 //cout<<"#################### Next event ##################"<<endl;
321 if(fRunHist) RunESDhist();
322 if(fRunTree) RunESDtree();
327 if(fRunHist) RunAODhist();
328 if(fRunTree) RunAODtree();
332 //_____________________________________________________________________________
333 void AliAnalysisTaskUpcPsi2s::RunAODtrig()
337 AliAODEvent *aod = (AliAODEvent*) InputEvent();
340 fRunNum = aod ->GetRunNumber();
342 TString trigger = aod->GetFiredTriggerClasses();
344 if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
346 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
347 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
349 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
350 if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
352 //MB, Central and SemiCentral triggers
353 AliCentrality *centrality = aod->GetCentrality();
354 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
356 Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
357 //Double_t percentile = centrality->GetCentralityPercentile("V0M");
359 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
361 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0) fHistCentralTriggersPerRun->Fill(fRunNum);
363 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
365 PostData(3, fListTrig);
368 //_____________________________________________________________________________
369 void AliAnalysisTaskUpcPsi2s::RunAODhist()
372 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
374 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
375 Double_t muonMass = partMuon->Mass();
377 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
378 Double_t electronMass = partElectron->Mass();
380 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
381 Double_t pionMass = partPion->Mass();
384 AliAODEvent *aod = (AliAODEvent*) InputEvent();
387 fHistNeventsJPsi->Fill(1);
388 fHistNeventsPsi2s->Fill(1);
391 TString trigger = aod->GetFiredTriggerClasses();
393 if( !trigger.Contains("CCUP") ) return;
395 fHistNeventsJPsi->Fill(2);
396 fHistNeventsPsi2s->Fill(2);
399 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
400 fVtxContrib = fAODVertex->GetNContributors();
401 if(fVtxContrib < 2) return;
403 fHistNeventsJPsi->Fill(3);
404 fHistNeventsPsi2s->Fill(3);
407 AliAODVZERO *fV0data = aod ->GetVZEROData();
408 AliAODZDC *fZDCdata = aod->GetZDCData();
410 fV0Adecision = fV0data->GetV0ADecision();
411 fV0Cdecision = fV0data->GetV0CDecision();
412 if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
414 fHistNeventsJPsi->Fill(4);
415 fHistNeventsPsi2s->Fill(4);
417 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
418 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
420 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
422 fHistNeventsJPsi->Fill(5);
423 fHistNeventsPsi2s->Fill(5);
426 Int_t nGoodTracks = 0;
427 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
429 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
430 Short_t qLepton[4], qPion[4];
431 UInt_t nLepton=0, nPion=0, nHighPt=0;
432 Double_t fRecTPCsignal[5];
433 Int_t mass[3]={-1,-1,-1};
437 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
438 AliAODTrack *trk = aod->GetTrack(itr);
440 if(!(trk->TestFilterBit(1<<0))) continue;
442 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
443 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
444 if(trk->GetTPCNcls() < 50)continue;
445 if(trk->Chi2perNDF() > 4)continue;
446 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
447 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
448 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
451 if(TMath::Abs(dca[1]) > 2) continue;
453 TrackIndex[nGoodTracks] = itr;
456 if(nGoodTracks > 4) break;
459 nLepton=0; nPion=0; nHighPt=0;
460 mass[0]= -1; mass[1]= -1, mass[2]= -1;
462 if(nGoodTracks == 4){
463 fHistNeventsPsi2s->Fill(6);
464 for(Int_t i=0; i<4; i++){
465 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
468 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
469 qLepton[nLepton] = trk->Charge();
470 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
471 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
474 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
475 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
481 qPion[nPion] = trk->Charge();
482 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
486 if(nLepton > 2 || nPion > 2) break;
488 if((nLepton == 2) && (nPion == 2)){
489 fHistNeventsPsi2s->Fill(7);
490 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(8);
491 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(9);
492 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
493 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
494 fHistNeventsPsi2s->Fill(11);
495 if(mass[0] == mass[1]) {
496 fHistNeventsPsi2s->Fill(12);
497 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
498 vDilepton = vLepton[0]+vLepton[1];
499 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
500 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
501 if(mass[0] == 0) fHistNeventsPsi2s->Fill(13);
502 if(mass[0] == 1) fHistNeventsPsi2s->Fill(14);
510 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
511 AliAODTrack *trk = aod->GetTrack(itr);
513 if(!(trk->TestFilterBit(1<<0))) continue;
515 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
516 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
517 if(trk->GetTPCNcls() < 70)continue;
518 if(trk->Chi2perNDF() > 4)continue;
519 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
520 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
521 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
522 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
524 if(TMath::Abs(dca[1]) > 2) continue;
525 if(TMath::Abs(dca[0]) > 0.2) continue;
527 TrackIndex[nGoodTracks] = itr;
530 if(nGoodTracks > 2) break;
533 nLepton=0; nPion=0; nHighPt=0;
534 mass[0]= -1; mass[1]= -1, mass[2]= -1;
536 if(nGoodTracks == 2){
537 fHistNeventsJPsi->Fill(6);
538 for(Int_t i=0; i<2; i++){
539 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
540 if(trk->Pt() > 1) nHighPt++;
541 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
542 qLepton[nLepton] = trk->Charge();
543 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
544 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
547 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
548 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
554 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(7);
555 if(qLepton[0]*qLepton[1] < 0){
556 fHistNeventsJPsi->Fill(8);
558 fHistNeventsJPsi->Fill(9);
559 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
560 if(nHighPt == 2) fHistNeventsJPsi->Fill(10);
561 if(mass[0] == mass[1] && mass[0] != -1) {
562 fHistNeventsJPsi->Fill(11);
563 vCandidate = vLepton[0]+vLepton[1];
564 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
566 fHistDiMuonMass->Fill(vCandidate.M());
567 fHistNeventsJPsi->Fill(12);
570 fHistDiElectronMass->Fill(vCandidate.M());
571 fHistNeventsJPsi->Fill(13);
580 PostData(4, fListHist);
584 //_____________________________________________________________________________
585 void AliAnalysisTaskUpcPsi2s::RunAODtree()
588 AliAODEvent *aod = (AliAODEvent*) InputEvent();
592 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
593 fDataFilnam->Clear();
594 fDataFilnam->SetString(filnam);
595 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
596 fRunNum = aod ->GetRunNumber();
599 TString trigger = aod->GetFiredTriggerClasses();
601 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
602 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
603 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
605 Bool_t isTriggered = kFALSE;
606 for(Int_t i=0; i<ntrg; i++) {
607 if( fTrigger[i] ) isTriggered = kTRUE;
609 if( !isTriggered ) return;
612 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
613 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
615 //Event identification
616 fPerNum = aod ->GetPeriodNumber();
617 fOrbNum = aod ->GetOrbitNumber();
618 fBCrossNum = aod ->GetBunchCrossNumber();
621 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
622 fVtxContrib = fAODVertex->GetNContributors();
623 fVtxPos[0] = fAODVertex->GetX();
624 fVtxPos[1] = fAODVertex->GetY();
625 fVtxPos[2] = fAODVertex->GetZ();
627 fAODVertex->GetCovarianceMatrix(CovMatx);
628 fVtxErr[0] = CovMatx[0];
629 fVtxErr[1] = CovMatx[1];
630 fVtxErr[2] = CovMatx[2];
631 fVtxChi2 = fAODVertex->GetChi2();
632 fVtxNDF = fAODVertex->GetNDF();
635 fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
638 AliAODVZERO *fV0data = aod ->GetVZEROData();
639 AliAODZDC *fZDCdata = aod->GetZDCData();
641 fV0Adecision = fV0data->GetV0ADecision();
642 fV0Cdecision = fV0data->GetV0CDecision();
643 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
644 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
648 //Track loop - loose cuts
649 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
650 AliAODTrack *trk = aod->GetTrack(itr);
652 if(!(trk->TestFilterBit(1<<0))) continue;
654 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
655 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
656 if(trk->GetTPCNcls() < 20)continue;
658 }//Track loop -loose cuts
661 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
664 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
665 AliAODTrack *trk = aod->GetTrack(itr);
667 if(!(trk->TestFilterBit(1<<0))) continue;
669 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
670 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
671 if(trk->GetTPCNcls() < 70)continue;
672 if(trk->Chi2perNDF() > 4)continue;
673 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
674 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
675 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
676 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
678 if(TMath::Abs(dca[1]) > 2) continue;
679 if(TMath::Abs(dca[0]) > 0.2) continue;
681 TrackIndex[nGoodTracks] = itr;
684 if(nGoodTracks > 2) break;
687 if(nGoodTracks == 2){
689 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
690 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
691 Double_t muonMass = partMuon->Mass();
692 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
693 Double_t electronMass = partElectron->Mass();
694 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
695 Double_t pionMass = partPion->Mass();
699 Double_t KFmass = pionMass;
700 Double_t fRecTPCsignal;
701 AliKFParticle *KFpart[2];
702 AliKFVertex *KFvtx = new AliKFVertex();
703 KFvtx->SetField(aod->GetMagneticField());
705 for(Int_t i=0; i<2; i++){
706 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
708 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
709 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
710 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
713 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
714 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
716 trk->GetPosition(KFpar);
717 trk->PxPyPz(KFpar+3);
718 trk->GetCovMatrix(KFcov);
721 fRecTPCsignal = trk->GetTPCsignal();
722 if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
723 if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
725 else KFmass = pionMass;
727 KFpart[i] = new AliKFParticle();
728 KFpart[i]->SetField(aod->GetMagneticField());
729 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
730 KFvtx->AddDaughter(*KFpart[i]);
733 Double_t pos[3]={0,0,0};
734 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
735 parTrk->CopyFromVTrack((AliVTrack*) trk);
736 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
738 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
739 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
743 fKfVtxPos[0]= KFvtx->GetX();
744 fKfVtxPos[1]= KFvtx->GetY();
745 fKfVtxPos[2]= KFvtx->GetZ();
746 for(UInt_t i=0; i<2; i++)delete KFpart[i];
754 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
755 AliAODTrack *trk = aod->GetTrack(itr);
757 if(!(trk->TestFilterBit(1<<0))) continue;
759 if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
760 if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
761 if(trk->GetTPCNcls() < 50)continue;
762 if(trk->Chi2perNDF() > 4)continue;
763 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
764 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
765 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
767 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
768 if(TMath::Abs(dca[1]) > 2) continue;
770 TrackIndex[nGoodTracks] = itr;
773 if(nGoodTracks > 4) break;
777 if(nGoodTracks == 4){
779 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
780 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
781 Double_t muonMass = partMuon->Mass();
782 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
783 Double_t electronMass = partElectron->Mass();
784 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
785 Double_t pionMass = partPion->Mass();
789 Double_t KFmass = pionMass;
790 Double_t fRecTPCsignal;
791 AliKFParticle *KFpart[4];
792 AliKFVertex *KFvtx = new AliKFVertex();
793 KFvtx->SetField(aod->GetMagneticField());
795 for(Int_t i=0; i<4; i++){
796 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
798 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
799 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
800 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
803 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
804 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
806 trk->GetPosition(KFpar);
807 trk->PxPyPz(KFpar+3);
808 trk->GetCovMatrix(KFcov);
811 fRecTPCsignal = trk->GetTPCsignal();
812 if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
813 if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
815 else KFmass = pionMass;
817 KFpart[i] = new AliKFParticle();
818 KFpart[i]->SetField(aod->GetMagneticField());
819 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
820 KFvtx->AddDaughter(*KFpart[i]);
822 Double_t pos[3]={0,0,0};
823 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
824 parTrk->CopyFromVTrack((AliVTrack*) trk);
825 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
827 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
828 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
832 fKfVtxPos[0]= KFvtx->GetX();
833 fKfVtxPos[1]= KFvtx->GetY();
834 fKfVtxPos[2]= KFvtx->GetZ();
835 for(UInt_t i=0; i<4; i++)delete KFpart[i];
840 PostData(1, fJPsiTree);
841 PostData(2, fPsi2sTree);
845 //_____________________________________________________________________________
846 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
850 AliESDEvent *esd = (AliESDEvent*) InputEvent();
853 fRunNum = esd ->GetRunNumber();
855 TString trigger = esd->GetFiredTriggerClasses();
857 if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
859 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
860 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
862 if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
864 //MB, Central and SemiCentral triggers
865 AliCentrality *centrality = esd->GetCentrality();
866 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
868 //Double_t percentile = centrality->GetCentralityPercentile("V0M");
869 Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
871 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
873 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0) fHistCentralTriggersPerRun->Fill(fRunNum);
875 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
878 PostData(3, fListTrig);
881 //_____________________________________________________________________________
882 void AliAnalysisTaskUpcPsi2s::RunESDhist()
885 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
887 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
888 Double_t muonMass = partMuon->Mass();
890 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
891 Double_t electronMass = partElectron->Mass();
893 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
894 Double_t pionMass = partPion->Mass();
897 AliESDEvent *esd = (AliESDEvent*) InputEvent();
900 fHistNeventsJPsi->Fill(1);
901 fHistNeventsPsi2s->Fill(1);
904 TString trigger = esd->GetFiredTriggerClasses();
906 if( !trigger.Contains("CCUP") ) return;
908 fHistNeventsJPsi->Fill(2);
909 fHistNeventsPsi2s->Fill(2);
912 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
913 fVtxContrib = fESDVertex->GetNContributors();
914 if(fVtxContrib < 2) return;
916 fHistNeventsJPsi->Fill(3);
917 fHistNeventsPsi2s->Fill(3);
920 AliESDVZERO *fV0data = esd->GetVZEROData();
921 AliESDZDC *fZDCdata = esd->GetESDZDC();
923 fV0Adecision = fV0data->GetV0ADecision();
924 fV0Cdecision = fV0data->GetV0CDecision();
925 if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
927 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
928 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
929 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
931 fHistNeventsJPsi->Fill(4);
932 fHistNeventsPsi2s->Fill(4);
936 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
938 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
939 Short_t qLepton[4], qPion[4];
940 UInt_t nLepton=0, nPion=0, nHighPt=0;
941 Double_t fRecTPCsignal[5];
942 Int_t mass[3]={-1,-1,-1};
945 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
946 AliESDtrack *trk = esd->GetTrack(itr);
949 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
950 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
951 if(trk->GetTPCNcls() < 70)continue;
952 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
953 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
954 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
955 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
956 trk->GetImpactParameters(dca[0],dca[1]);
957 if(TMath::Abs(dca[1]) > 2) continue;
958 if(TMath::Abs(dca[1]) > 0.2) continue;
960 TrackIndex[nGoodTracks] = itr;
962 if(nGoodTracks > 2) break;
966 if(nGoodTracks == 2){
967 fHistNeventsJPsi->Fill(5);
968 for(Int_t i=0; i<2; i++){
969 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
970 if(trk->Pt() > 1) nHighPt++;
971 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
972 qLepton[nLepton] = trk->Charge();
973 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
974 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
977 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
978 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
984 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
985 if(qLepton[0]*qLepton[1] < 0){
986 fHistNeventsJPsi->Fill(7);
988 fHistNeventsJPsi->Fill(8);
989 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
990 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
991 if(mass[0] == mass[1] && mass[0] != -1) {
992 fHistNeventsJPsi->Fill(10);
993 vCandidate = vLepton[0]+vLepton[1];
994 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
996 fHistDiMuonMass->Fill(vCandidate.M());
997 fHistNeventsJPsi->Fill(11);
1000 fHistDiElectronMass->Fill(vCandidate.M());
1001 fHistNeventsJPsi->Fill(12);
1008 nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1009 mass[0]= -1; mass[1]= -1, mass[2]= -1;
1012 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1013 AliESDtrack *trk = esd->GetTrack(itr);
1014 if( !trk ) continue;
1016 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1017 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1018 if(trk->GetTPCNcls() < 50)continue;
1019 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1020 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1021 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1022 trk->GetImpactParameters(dca[0],dca[1]);
1023 if(TMath::Abs(dca[1]) > 2) continue;
1025 TrackIndex[nGoodTracks] = itr;
1027 if(nGoodTracks > 4) break;
1030 if(nGoodTracks == 4){
1031 fHistNeventsPsi2s->Fill(5);
1032 for(Int_t i=0; i<4; i++){
1033 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1036 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1037 qLepton[nLepton] = trk->Charge();
1038 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1039 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1042 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1043 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1049 qPion[nPion] = trk->Charge();
1050 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1054 if(nLepton > 2 || nPion > 2) break;
1056 if((nLepton == 2) && (nPion == 2)){
1057 fHistNeventsPsi2s->Fill(6);
1058 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1059 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1060 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1061 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1062 fHistNeventsPsi2s->Fill(10);
1063 if(mass[0] == mass[1]) {
1064 fHistNeventsPsi2s->Fill(11);
1065 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1066 vDilepton = vLepton[0]+vLepton[1];
1067 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1068 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1069 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1070 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1076 PostData(4, fListHist);
1080 //_____________________________________________________________________________
1081 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1085 AliESDEvent *esd = (AliESDEvent*) InputEvent();
1089 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1090 fDataFilnam->Clear();
1091 fDataFilnam->SetString(filnam);
1092 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1093 fRunNum = esd->GetRunNumber();
1096 TString trigger = esd->GetFiredTriggerClasses();
1098 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1099 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
1100 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1102 Bool_t isTriggered = kFALSE;
1103 for(Int_t i=0; i<ntrg; i++) {
1104 if( fTrigger[i] ) isTriggered = kTRUE;
1106 if( !isTriggered ) return;
1109 fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1110 fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1112 //TOF trigger info (0OMU)
1113 fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1114 fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1116 //Event identification
1117 fPerNum = esd->GetPeriodNumber();
1118 fOrbNum = esd->GetOrbitNumber();
1119 fBCrossNum = esd->GetBunchCrossNumber();
1122 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1123 fVtxContrib = fESDVertex->GetNContributors();
1124 fVtxPos[0] = fESDVertex->GetX();
1125 fVtxPos[1] = fESDVertex->GetY();
1126 fVtxPos[2] = fESDVertex->GetZ();
1127 Double_t CovMatx[6];
1128 fESDVertex->GetCovarianceMatrix(CovMatx);
1129 fVtxErr[0] = CovMatx[0];
1130 fVtxErr[1] = CovMatx[1];
1131 fVtxErr[2] = CovMatx[2];
1132 fVtxChi2 = fESDVertex->GetChi2();
1133 fVtxNDF = fESDVertex->GetNDF();
1136 fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1139 AliESDVZERO *fV0data = esd->GetVZEROData();
1140 AliESDZDC *fZDCdata = esd->GetESDZDC();
1142 fV0Adecision = fV0data->GetV0ADecision();
1143 fV0Cdecision = fV0data->GetV0CDecision();
1144 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1145 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1149 //Track loop - loose cuts
1150 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1151 AliESDtrack *trk = esd->GetTrack(itr);
1152 if( !trk ) continue;
1154 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1155 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1156 if(trk->GetTPCNcls() < 20)continue;
1158 }//Track loop -loose cuts
1160 Int_t nGoodTracks=0;
1161 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1164 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1165 AliESDtrack *trk = esd->GetTrack(itr);
1166 if( !trk ) continue;
1168 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1169 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1170 if(trk->GetTPCNcls() < 70)continue;
1171 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1172 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1173 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1174 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1175 trk->GetImpactParameters(dca[0],dca[1]);
1176 if(TMath::Abs(dca[1]) > 2) continue;
1177 if(TMath::Abs(dca[1]) > 0.2) continue;
1179 TrackIndex[nGoodTracks] = itr;
1181 if(nGoodTracks > 2) break;
1184 if(nGoodTracks == 2){
1185 for(Int_t i=0; i<2; i++){
1186 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1188 AliExternalTrackParam cParam;
1189 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1191 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
1193 Double_t pos[3]={0,0,0};
1194 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1196 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1197 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1205 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1206 AliESDtrack *trk = esd->GetTrack(itr);
1207 if( !trk ) continue;
1209 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1210 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1211 if(trk->GetTPCNcls() < 50)continue;
1212 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1213 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1214 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1215 trk->GetImpactParameters(dca[0],dca[1]);
1216 if(TMath::Abs(dca[1]) > 2) continue;
1218 TrackIndex[nGoodTracks] = itr;
1220 if(nGoodTracks > 4) break;
1223 if(nGoodTracks == 4){
1224 for(Int_t i=0; i<4; i++){
1225 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1227 AliExternalTrackParam cParam;
1228 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1230 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1232 Double_t pos[3]={0,0,0};
1233 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1235 fTOFphi[i] = TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1236 if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1239 fPsi2sTree ->Fill();
1242 PostData(1, fJPsiTree);
1243 PostData(2, fPsi2sTree);
1247 //_____________________________________________________________________________
1248 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
1251 cout<<"Analysis complete."<<endl;