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"
53 #include "AliAnalysisTaskUpcPsi2s.h"
55 ClassImp(AliAnalysisTaskUpcPsi2s);
60 //trees for UPC analysis,
61 // michal.broz@cern.ch
63 //_____________________________________________________________________________
64 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s()
65 : AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fJPsiTree(0),fPsi2sTree(0),
66 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
67 fVtxContrib(0),fVtxPosX(0),fVtxPosY(0),fVtxPosZ(0),fVtxErrX(0),fVtxErrY(0),fVtxErrZ(0),fVtxChi2(0),fVtxNDF(0),fVtxType(0),
68 fBCrossNum(0),fNtracklets(0),
69 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
70 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
71 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
72 fListTrig(0),fHistUpcTriggersPerRun(0),fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0),
73 fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
74 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
75 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
81 }//AliAnalysisTaskUpcPsi2s
84 //_____________________________________________________________________________
85 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
86 : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fJPsiTree(0),fPsi2sTree(0),
87 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
88 fVtxContrib(0),fVtxPosX(0),fVtxPosY(0),fVtxPosZ(0),fVtxErrX(0),fVtxErrY(0),fVtxErrZ(0),fVtxChi2(0),fVtxNDF(0),fVtxType(0),
89 fBCrossNum(0),fNtracklets(0),
90 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
91 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
92 fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),
93 fListTrig(0),fHistUpcTriggersPerRun(0),fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0),
94 fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
95 fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
96 fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
101 if( strstr(name,"ESD") ) fType = 0;
102 if( strstr(name,"AOD") ) fType = 1;
106 DefineOutput(1, TTree::Class());
107 DefineOutput(2, TTree::Class());
108 DefineOutput(3, TList::Class());
109 DefineOutput(4, TList::Class());
111 }//AliAnalysisTaskUpcPsi2s
113 //_____________________________________________________________________________
114 void AliAnalysisTaskUpcPsi2s::Init()
117 for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
121 //_____________________________________________________________________________
122 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
142 }//~AliAnalysisTaskUpcPsi2s
145 //_____________________________________________________________________________
146 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
149 fDataFilnam = new TObjString();
150 fDataFilnam->SetString("");
153 fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
154 fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
155 fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
156 fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
158 //output tree with JPsi candidate events
159 fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
160 fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
161 fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
162 fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
164 fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
165 fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
166 fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
167 fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
168 fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
169 fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
171 fJPsiTree ->Branch("fVtxPosX", &fVtxPosX, "fVtxPosX/D");
172 fJPsiTree ->Branch("fVtxPosY", &fVtxPosY, "fVtxPosY/D");
173 fJPsiTree ->Branch("fVtxPosZ", &fVtxPosZ, "fVtxPosZ/D");
174 fJPsiTree ->Branch("fVtxErrX", &fVtxErrX, "fVtxErrX/D");
175 fJPsiTree ->Branch("fVtxErrY", &fVtxErrY, "fVtxErrY/D");
176 fJPsiTree ->Branch("fVtxErrZ", &fVtxErrZ, "fVtxErrZ/D");
177 fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
178 fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
179 fJPsiTree ->Branch("fVtxType", &fVtxType, "fVtxType/B");
181 fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
182 fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
183 fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
184 fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
185 fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
186 fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
187 fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
189 fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
192 fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
195 //output tree with Psi2s candidate events
196 fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
197 fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
198 fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
199 fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
201 fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
202 fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
203 fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
204 fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
205 fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
206 fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
208 fPsi2sTree ->Branch("fVtxPosX", &fVtxPosX, "fVtxPosX/D");
209 fPsi2sTree ->Branch("fVtxPosY", &fVtxPosY, "fVtxPosY/D");
210 fPsi2sTree ->Branch("fVtxPosZ", &fVtxPosZ, "fVtxPosZ/D");
211 fPsi2sTree ->Branch("fVtxErrX", &fVtxErrX, "fVtxErrX/D");
212 fPsi2sTree ->Branch("fVtxErrY", &fVtxErrY, "fVtxErrY/D");
213 fPsi2sTree ->Branch("fVtxErrZ", &fVtxErrZ, "fVtxErrZ/D");
214 fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
215 fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
216 fPsi2sTree ->Branch("fVtxType", &fVtxType, "fVtxType/B");
218 fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
219 fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
220 fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
221 fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
222 fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
223 fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
224 fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
226 fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
229 fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
232 fListTrig = new TList();
233 fListTrig ->SetOwner();
235 fHistUpcTriggersPerRun = new TH1D("fHistUpcTriggersPerRun", "fHistUpcTriggersPerRun", 3000, 167000.5, 170000.5);
236 fListTrig->Add(fHistUpcTriggersPerRun);
238 fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 3000, 167000.5, 170000.5);
239 fListTrig->Add(fHistZedTriggersPerRun);
241 fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 3000, 167000.5, 170000.5);
242 fListTrig->Add(fHistCvlnTriggersPerRun);
244 fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 3000, 167000.5, 170000.5);
245 fListTrig->Add(fHistMBTriggersPerRun);
247 fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 3000, 167000.5, 170000.5);
248 fListTrig->Add(fHistCentralTriggersPerRun);
250 fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 3000, 167000.5, 170000.5);
251 fListTrig->Add(fHistSemiCentralTriggersPerRun);
253 fListHist = new TList();
254 fListHist ->SetOwner();
256 TString CutNameJPsi[12] = {"Analyzed","Triggered","Vertex cut","V0 decision","Two good tracks",
257 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
258 fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",12,0.5,12.5);
259 for (Int_t i = 0; i<12; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
260 fListHist->Add(fHistNeventsJPsi);
262 fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
263 fListHist->Add(fHistTPCsignalJPsi);
265 fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
266 fListHist->Add(fHistDiLeptonPtJPsi);
268 fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
269 fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
270 fListHist->Add(fHistDiElectronMass);
272 fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
273 fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
274 fListHist->Add(fHistDiMuonMass);
276 TString CutNamePsi2s[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Four good tracks",
277 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
279 fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",13,0.5,13.5);
280 for (Int_t i = 0; i<13; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
281 fListHist->Add(fHistNeventsPsi2s);
283 fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
284 fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
285 fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
286 fListHist->Add(fHistPsi2sMassVsPt);
288 fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
289 fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
290 fListHist->Add(fHistPsi2sMassCoherent);
292 PostData(1, fJPsiTree);
293 PostData(2, fPsi2sTree);
294 PostData(3, fListTrig);
295 PostData(4, fListHist);
297 }//UserCreateOutputObjects
299 //_____________________________________________________________________________
300 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
303 //cout<<"#################### Next event ##################"<<endl;
307 if(fRunHist) RunESDhist();
308 if(fRunTree) RunESDtree();
313 if(fRunHist) RunAODhist();
314 if(fRunTree) RunAODtree();
318 //_____________________________________________________________________________
319 void AliAnalysisTaskUpcPsi2s::RunAODtrig()
323 AliAODEvent *aod = (AliAODEvent*) InputEvent();
326 fRunNum = aod ->GetRunNumber();
328 TString trigger = aod->GetFiredTriggerClasses();
330 if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
332 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
333 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
335 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
336 if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
338 //MB, Central and SemiCentral triggers
339 AliCentrality *centrality = aod->GetCentrality();
340 UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
342 Double_t percentile = centrality->GetCentralityPercentile("V0M");
344 if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<80 && percentile>0) fHistMBTriggersPerRun->Fill(fRunNum);
346 if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<6 && percentile>0) fHistCentralTriggersPerRun->Fill(fRunNum);
348 if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<50 && percentile>15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
350 PostData(3, fListTrig);
353 //_____________________________________________________________________________
354 void AliAnalysisTaskUpcPsi2s::RunAODhist()
357 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
359 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
360 Double_t muonMass = partMuon->Mass();
362 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
363 Double_t electronMass = partElectron->Mass();
365 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
366 Double_t pionMass = partPion->Mass();
369 AliAODEvent *aod = (AliAODEvent*) InputEvent();
372 fHistNeventsJPsi->Fill(1);
373 fHistNeventsPsi2s->Fill(1);
376 TString trigger = aod->GetFiredTriggerClasses();
378 if( !trigger.Contains("CCUP4-B") ) return;
380 fHistNeventsJPsi->Fill(2);
381 fHistNeventsPsi2s->Fill(2);
384 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
385 fVtxContrib = fAODVertex->GetNContributors();
386 if(fVtxContrib < 2) return;
388 fHistNeventsJPsi->Fill(3);
389 fHistNeventsPsi2s->Fill(3);
392 AliAODVZERO *fV0data = aod ->GetVZEROData();
393 AliAODZDC *fZDCdata = aod->GetZDCData();
395 fV0Adecision = fV0data->GetV0ADecision();
396 fV0Cdecision = fV0data->GetV0CDecision();
397 if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
399 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
400 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
402 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
404 fHistNeventsJPsi->Fill(4);
405 fHistNeventsPsi2s->Fill(4);
409 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
411 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
412 Short_t qLepton[4], qPion[4];
413 UInt_t nLepton=0, nPion=0, nHighPt=0;
414 Double_t jRecTPCsignal[5];
415 Int_t mass[3]={-1,-1,-1};
418 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
419 AliAODTrack *trk = aod->GetTrack(itr);
422 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
423 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
424 if(trk->GetTPCNcls() < 70)continue;
425 if(trk->Chi2perNDF() > 4)continue;
426 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
427 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
428 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
429 if(TMath::Abs(dca[1]) > 2) continue;
430 if(TMath::Abs(dca[0]) > 0.2) continue;
432 TrackIndex[nGoodTracks] = itr;
435 if(nGoodTracks > 2) break;
438 if(nGoodTracks == 2){
439 fHistNeventsJPsi->Fill(5);
440 for(Int_t i=0; i<2; i++){
441 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
442 if(trk->Pt() > 1) nHighPt++;
443 jRecTPCsignal[nLepton] = trk->GetTPCsignal();
444 qLepton[nLepton] = trk->Charge();
445 if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
446 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
449 if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
450 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
456 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
457 if(qLepton[0]*qLepton[1] < 0){
458 fHistNeventsJPsi->Fill(7);
460 fHistNeventsJPsi->Fill(8);
461 fHistTPCsignalJPsi->Fill(jRecTPCsignal[0],jRecTPCsignal[1]);
462 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
463 if(mass[0] == mass[1] && mass[0] != -1) {
464 fHistNeventsJPsi->Fill(10);
465 vCandidate = vLepton[0]+vLepton[1];
466 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
468 fHistDiMuonMass->Fill(vCandidate.M());
469 fHistNeventsJPsi->Fill(11);
472 fHistDiElectronMass->Fill(vCandidate.M());
473 fHistNeventsJPsi->Fill(12);
483 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
484 AliAODTrack *trk = aod->GetTrack(itr);
487 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
488 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
489 if(trk->GetTPCNcls() < 50)continue;
490 if(trk->Chi2perNDF() > 4)continue;
491 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
492 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
493 if(TMath::Abs(dca[1]) > 2) continue;
495 TrackIndex[nGoodTracks] = itr;
498 if(nGoodTracks > 4) break;
501 nLepton=0; nPion=0; nHighPt=0;
502 mass[0]= -1; mass[1]= -1, mass[2]= -1;
504 if(nGoodTracks == 4){
505 fHistNeventsPsi2s->Fill(5);
506 for(Int_t i=0; i<4; i++){
507 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
510 jRecTPCsignal[nLepton] = trk->GetTPCsignal();
511 qLepton[nLepton] = trk->Charge();
512 if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
513 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
516 if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
517 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
523 qPion[nPion] = trk->Charge();
524 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
528 if(nLepton > 2 || nPion > 2) break;
530 if((nLepton == 2) && (nPion == 2)){
531 fHistNeventsPsi2s->Fill(6);
532 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
533 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
534 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
535 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
536 fHistNeventsPsi2s->Fill(10);
537 if(mass[0] == mass[1]) {
538 fHistNeventsPsi2s->Fill(11);
539 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
540 vDilepton = vLepton[0]+vLepton[1];
541 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
542 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
543 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
544 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
550 PostData(4, fListHist);
554 //_____________________________________________________________________________
555 void AliAnalysisTaskUpcPsi2s::RunAODtree()
558 AliAODEvent *aod = (AliAODEvent*) InputEvent();
562 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
563 fDataFilnam->Clear();
564 fDataFilnam->SetString(filnam);
565 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
566 fRunNum = aod ->GetRunNumber();
569 TString trigger = aod->GetFiredTriggerClasses();
571 if( !trigger.Contains("CCUP4-B") ) return;
574 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
575 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
577 //Event identification
578 fPerNum = aod ->GetPeriodNumber();
579 fOrbNum = aod ->GetOrbitNumber();
580 fBCrossNum = aod ->GetBunchCrossNumber();
583 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
584 fVtxContrib = fAODVertex->GetNContributors();
585 fVtxPosX = fAODVertex->GetX();
586 fVtxPosY = fAODVertex->GetY();
587 fVtxPosZ = fAODVertex->GetZ();
589 fAODVertex->GetCovarianceMatrix(CovMatx);
590 fVtxErrX = CovMatx[0];
591 fVtxErrY = CovMatx[1];
592 fVtxErrZ = CovMatx[2];
593 fVtxChi2 = fAODVertex->GetChi2();
594 fVtxNDF = fAODVertex->GetNDF();
595 fVtxType = fAODVertex->GetType();
598 fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
601 AliAODVZERO *fV0data = aod ->GetVZEROData();
602 AliAODZDC *fZDCdata = aod->GetZDCData();
604 fV0Adecision = fV0data->GetV0ADecision();
605 fV0Cdecision = fV0data->GetV0CDecision();
606 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
607 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
610 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
613 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
614 AliAODTrack *trk = aod->GetTrack(itr);
617 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
618 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
619 if(trk->GetTPCNcls() < 70)continue;
620 if(trk->Chi2perNDF() > 4)continue;
621 if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
622 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
623 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
624 if(TMath::Abs(dca[1]) > 2) continue;
625 if(TMath::Abs(dca[0]) > 0.2) continue;
627 TrackIndex[nGoodTracks] = itr;
630 if(nGoodTracks > 2) break;
633 if(nGoodTracks == 2){
634 for(Int_t i=0; i<2; i++){
635 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
637 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
638 trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
640 trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
641 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
645 PostData(1, fJPsiTree);
650 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
651 AliAODTrack *trk = aod->GetTrack(itr);
654 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
655 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
656 if(trk->GetTPCNcls() < 50)continue;
657 if(trk->Chi2perNDF() > 4)continue;
658 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
659 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
660 if(TMath::Abs(dca[1]) > 2) continue;
662 TrackIndex[nGoodTracks] = itr;
665 if(nGoodTracks > 4) break;
669 if(nGoodTracks == 4){
670 for(Int_t i=0; i<4; i++){
671 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
673 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
674 trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov);
676 trk->SetDCA(dca[0],dca[1]); //to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
677 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
681 PostData(2, fPsi2sTree);
686 //_____________________________________________________________________________
687 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
691 AliESDEvent *esd = (AliESDEvent*) InputEvent();
694 fRunNum = esd ->GetRunNumber();
696 TString trigger = esd->GetFiredTriggerClasses();
698 if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
700 if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
701 if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
703 if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
705 PostData(3, fListTrig);
708 //_____________________________________________________________________________
709 void AliAnalysisTaskUpcPsi2s::RunESDhist()
712 TDatabasePDG *pdgdat = TDatabasePDG::Instance();
714 TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
715 Double_t muonMass = partMuon->Mass();
717 TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
718 Double_t electronMass = partElectron->Mass();
720 TParticlePDG *partPion = pdgdat->GetParticle( 211 );
721 Double_t pionMass = partPion->Mass();
724 AliESDEvent *esd = (AliESDEvent*) InputEvent();
727 fHistNeventsJPsi->Fill(1);
728 fHistNeventsPsi2s->Fill(1);
731 TString trigger = esd->GetFiredTriggerClasses();
733 if( !trigger.Contains("CCUP4-B") ) return;
735 fHistNeventsJPsi->Fill(2);
736 fHistNeventsPsi2s->Fill(2);
739 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
740 fVtxContrib = fESDVertex->GetNContributors();
741 if(fVtxContrib < 2) return;
743 fHistNeventsJPsi->Fill(3);
744 fHistNeventsPsi2s->Fill(3);
747 AliESDVZERO *fV0data = esd->GetVZEROData();
748 AliESDZDC *fZDCdata = esd->GetESDZDC();
750 fV0Adecision = fV0data->GetV0ADecision();
751 fV0Cdecision = fV0data->GetV0CDecision();
752 if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
754 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
755 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
756 if( fZDCAenergy > 12000 || fZDCCenergy > 12000) return;
758 fHistNeventsJPsi->Fill(4);
759 fHistNeventsPsi2s->Fill(4);
763 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
765 TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
766 Short_t qLepton[4], qPion[4];
767 UInt_t nLepton=0, nPion=0, nHighPt=0;
768 Double_t jRecTPCsignal[5];
769 Int_t mass[3]={-1,-1,-1};
772 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
773 AliESDtrack *trk = esd->GetTrack(itr);
776 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
777 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
778 if(trk->GetTPCNcls() < 50)continue;
779 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
780 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
781 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
782 trk->GetImpactParameters(dca[0],dca[1]);
783 if(TMath::Abs(dca[1]) > 2) continue;
785 TrackIndex[nGoodTracks] = itr;
787 if(nGoodTracks > 4) break;
790 if(nGoodTracks == 2){
791 fHistNeventsJPsi->Fill(5);
792 for(Int_t i=0; i<2; i++){
793 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
794 if(trk->Pt() > 1) nHighPt++;
795 jRecTPCsignal[nLepton] = trk->GetTPCsignal();
796 qLepton[nLepton] = trk->Charge();
797 if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
798 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
801 if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
802 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
808 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
809 if(qLepton[0]*qLepton[1] < 0){
810 fHistNeventsJPsi->Fill(7);
812 fHistNeventsJPsi->Fill(8);
813 fHistTPCsignalJPsi->Fill(jRecTPCsignal[0],jRecTPCsignal[1]);
814 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
815 if(mass[0] == mass[1] && mass[0] != -1) {
816 fHistNeventsJPsi->Fill(10);
817 vCandidate = vLepton[0]+vLepton[1];
818 if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
820 fHistDiMuonMass->Fill(vCandidate.M());
821 fHistNeventsJPsi->Fill(11);
824 fHistDiElectronMass->Fill(vCandidate.M());
825 fHistNeventsJPsi->Fill(12);
832 nLepton=0; nPion=0; nHighPt=0;
833 mass[0]= -1; mass[1]= -1, mass[2]= -1;
835 if(nGoodTracks == 4){
836 fHistNeventsPsi2s->Fill(5);
837 for(Int_t i=0; i<4; i++){
838 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
841 jRecTPCsignal[nLepton] = trk->GetTPCsignal();
842 qLepton[nLepton] = trk->Charge();
843 if(jRecTPCsignal[nLepton] > 40 && jRecTPCsignal[nLepton] < 70){
844 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
847 if(jRecTPCsignal[nLepton] > 70 && jRecTPCsignal[nLepton] < 100){
848 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
854 qPion[nPion] = trk->Charge();
855 vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
859 if(nLepton > 2 || nPion > 2) break;
861 if((nLepton == 2) && (nPion == 2)){
862 fHistNeventsPsi2s->Fill(6);
863 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
864 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
865 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
866 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
867 fHistNeventsPsi2s->Fill(10);
868 if(mass[0] == mass[1]) {
869 fHistNeventsPsi2s->Fill(11);
870 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
871 vDilepton = vLepton[0]+vLepton[1];
872 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
873 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
874 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
875 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
881 PostData(4, fListHist);
885 //_____________________________________________________________________________
886 void AliAnalysisTaskUpcPsi2s::RunESDtree()
890 AliESDEvent *esd = (AliESDEvent*) InputEvent();
894 const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
895 fDataFilnam->Clear();
896 fDataFilnam->SetString(filnam);
897 fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
898 fRunNum = esd->GetRunNumber();
901 TString trigger = esd->GetFiredTriggerClasses();
903 if( !trigger.Contains("CCUP4-B") ) return;
906 fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
907 fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
909 //Event identification
910 fPerNum = esd->GetPeriodNumber();
911 fOrbNum = esd->GetOrbitNumber();
912 fBCrossNum = esd->GetBunchCrossNumber();
915 AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
916 fVtxContrib = fESDVertex->GetNContributors();
919 fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
922 AliESDVZERO *fV0data = esd->GetVZEROData();
923 AliESDZDC *fZDCdata = esd->GetESDZDC();
925 fV0Adecision = fV0data->GetV0ADecision();
926 fV0Cdecision = fV0data->GetV0CDecision();
927 fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
928 fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
931 Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
934 for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
935 AliESDtrack *trk = esd->GetTrack(itr);
938 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
939 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
940 if(trk->GetTPCNcls() < 50)continue;
941 if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
942 Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
943 if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
944 trk->GetImpactParameters(dca[0],dca[1]);
945 if(TMath::Abs(dca[1]) > 2) continue;
947 TrackIndex[nGoodTracks] = itr;
949 if(nGoodTracks > 4) break;
952 if(nGoodTracks == 2){
953 for(Int_t i=0; i<2; i++){
954 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
956 AliExternalTrackParam cParam;
957 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
959 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
963 PostData(1, fJPsiTree);
966 if(nGoodTracks == 4){
967 for(Int_t i=0; i<4; i++){
968 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
970 AliExternalTrackParam cParam;
971 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
973 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
977 PostData(2, fPsi2sTree);
982 //_____________________________________________________________________________
983 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
986 cout<<"Analysis complete."<<endl;