PWGUD/UPC/AliAnalysisTaskUpcPsi2s.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 // c++ headers
  17 #include <iostream>
  18 #include <string.h>
  19
  20 // root headers
  21 #include "TH1I.h"
  22 #include "TTree.h"
  23 #include "TClonesArray.h"
  24 #include "TParticle.h"
  25 #include "TObjString.h"
  26 #include "TFile.h"
  27 #include "TDatabasePDG.h"
  28 #include "TLorentzVector.h"
  29
  30 // aliroot headers
  31 #include "AliAnalysisManager.h"
  32 #include "AliInputEventHandler.h"
  33 #include "AliESDEvent.h"
  34 #include "AliAODEvent.h"
  35 #include "AliMCEvent.h"
  36 #include "AliAODVZERO.h"
  37 #include "AliAODZDC.h"
  38 #include "AliESDVZERO.h"
  39 #include "AliESDZDC.h"
  40 #include "AliPIDResponse.h"
  41 #include "AliAODTrack.h"
  42 #include "AliAODPid.h"
  43 #include "AliAODVertex.h"
  44 #include "AliESDVertex.h"
  45 #include "AliMultiplicity.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDMuonTrack.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliMCParticle.h"
  50 #include "AliCentrality.h"
  51 #include "AliKFVertex.h"
  52 #include "AliExternalTrackParam.h"
  53 #include "AliTriggerAnalysis.h"
  54
  55 // my headers
  56 #include "AliAnalysisTaskUpcPsi2s.h"
  57
  58 ClassImp(AliAnalysisTaskUpcPsi2s);
  59
  60 using std::cout;
  61 using std::endl;
  62
  63 //trees for UPC analysis,
  64 // michal.broz@cern.ch
  65
  66 //_____________________________________________________________________________
  67 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s()
  68   : AliAnalysisTaskSE(),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fRunSystematics(kFALSE),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
  69     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
  70     fTOFtrig1(0), fTOFtrig2(0),
  71     fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
  72     fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
  73     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
  74     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
  75     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
  76     fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),
  77     fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
  78     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),fHistDiLeptonMass(0),
  79     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
  80     fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0)
  81
  82 {
  83
  84 //Dummy constructor
  85
  86 }//AliAnalysisTaskUpcPsi2s
  87
  88
  89 //_____________________________________________________________________________
  90 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
  91   : AliAnalysisTaskSE(name),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fRunSystematics(kFALSE),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
  92     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
  93     fTOFtrig1(0), fTOFtrig2(0),
  94     fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
  95     fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
  96     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
  97     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
  98     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
  99     fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),
 100     fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
 101     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),fHistDiLeptonMass(0),
 102     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0),
 103     fListSystematics(0),fListJPsiLoose(0),fListJPsiTight(0)
 104
 105 {
 106
 107   // Constructor
 108   if( strstr(name,"ESD") ) fType = 0;
 109   if( strstr(name,"AOD") ) fType = 1;
 110
 111   Init();
 112
 113   DefineOutput(1, TTree::Class());
 114   DefineOutput(2, TTree::Class());
 115   DefineOutput(3, TList::Class());
 116   DefineOutput(4, TList::Class());
 117
 118 }//AliAnalysisTaskUpcPsi2s
 119
 120 //_____________________________________________________________________________
 121 void AliAnalysisTaskUpcPsi2s::Init()
 122 {
 123
 124   for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
 125   for(Int_t i=0; i<4; i++) {
 126         fTOFphi[i] = -666;
 127         fPIDMuon[i] = -666;
 128         fPIDElectron[i] = -666;
 129         fPIDPion[i] = -666;
 130         fTriggerInputsMC[i] = kFALSE;
 131         }
 132   for(Int_t i=0; i<3; i++){
 133         fVtxPos[i] = -666;
 134         fVtxErr[i] = -666;
 135         fKfVtxPos[i] = -666;
 136         }
 137
 138 }//Init
 139
 140 //_____________________________________________________________________________
 141 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
 142 {
 143   // Destructor
 144   if(fJPsiTree){
 145      delete fJPsiTree;
 146      fJPsiTree = 0x0;
 147   }
 148   if(fPsi2sTree){
 149      delete fPsi2sTree;
 150      fPsi2sTree = 0x0;
 151   }
 152   if(fListTrig){
 153      delete fListTrig;
 154      fListTrig = 0x0;
 155   }
 156   if(fListHist){
 157      delete fListHist;
 158      fListHist = 0x0;
 159   }
 160
 161 }//~AliAnalysisTaskUpcPsi2s
 162
 163
 164 //_____________________________________________________________________________
 165 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
 166 {
 167   //PID response
 168   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 169   AliInputEventHandler *inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 170   fPIDResponse = inputHandler->GetPIDResponse();
 171
 172   //input file
 173   fDataFilnam = new TObjString();
 174   fDataFilnam->SetString("");
 175
 176     //tracks
 177   fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
 178   fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
 179   fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
 180   fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
 181   fGenPart = new TClonesArray("TParticle", 1000);
 182
 183   //output tree with JPsi candidate events
 184   fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
 185   fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 186   fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 187   fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 188
 189   fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 190   fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 191   fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 192   fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 193   fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 194   fJPsiTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
 195   fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 196
 197   fJPsiTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
 198   fJPsiTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
 199   fJPsiTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
 200
 201   fJPsiTree ->Branch("fPIDMuon", &fPIDMuon[0], "fPIDMuon[4]/D");
 202   fJPsiTree ->Branch("fPIDElectron", &fPIDElectron[0], "fPIDElectron[4]/D");
 203   fJPsiTree ->Branch("fPIDPion", &fPIDPion[0], "fPIDPion[4]/D");
 204
 205   fJPsiTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
 206   fJPsiTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
 207   fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
 208   fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
 209
 210   fJPsiTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
 211
 212   fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 213   fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 214   fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 215   fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 216   fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
 217   fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 218   fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 219   if( fType == 0 ) {
 220     fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
 221   }
 222   if( fType == 1 ) {
 223     fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
 224   }
 225   if(isMC) {
 226     fJPsiTree ->Branch("fGenPart", &fGenPart);
 227     fJPsiTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
 228   }
 229
 230
 231  //output tree with Psi2s candidate events
 232   fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
 233   fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 234   fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 235   fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 236
 237   fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 238   fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 239   fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 240   fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 241   fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 242   fPsi2sTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
 243   fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 244
 245   fPsi2sTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
 246   fPsi2sTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
 247   fPsi2sTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
 248
 249   fPsi2sTree ->Branch("fPIDMuon", &fPIDMuon[0], "fPIDMuon[4]/D");
 250   fPsi2sTree ->Branch("fPIDElectron", &fPIDElectron[0], "fPIDElectron[4]/D");
 251   fPsi2sTree ->Branch("fPIDPion", &fPIDPion[0], "fPIDPion[4]/D");
 252
 253   fPsi2sTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
 254   fPsi2sTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
 255   fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
 256   fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
 257
 258   fPsi2sTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
 259
 260   fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 261   fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 262   fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 263   fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 264   fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
 265   fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 266   fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 267   if( fType == 0 ) {
 268     fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
 269   }
 270   if( fType == 1 ) {
 271     fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
 272   }
 273   if(isMC) {
 274     fPsi2sTree ->Branch("fGenPart", &fGenPart);
 275     fPsi2sTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
 276   }
 277
 278
 279   fListTrig = new TList();
 280   fListTrig ->SetOwner();
 281
 282   fHistCcup4TriggersPerRun = new TH1D("fHistCcup4TriggersPerRun", "fHistCcup4TriggersPerRun", 33000, 167000.5, 200000.5);
 283   fListTrig->Add(fHistCcup4TriggersPerRun);
 284
 285   fHistCcup7TriggersPerRun = new TH1D("fHistCcup7TriggersPerRun", "fHistCcup7TriggersPerRun", 33000, 167000.5, 200000.5);
 286   fListTrig->Add(fHistCcup7TriggersPerRun);
 287
 288   fHistCcup2TriggersPerRun = new TH1D("fHistCcup2TriggersPerRun", "fHistCcup2TriggersPerRun", 33000, 167000.5, 200000.5);
 289   fListTrig->Add(fHistCcup2TriggersPerRun);
 290
 291   fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 33000, 167000.5, 200000.5);
 292   fListTrig->Add(fHistZedTriggersPerRun);
 293
 294   fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 33000, 167000.5, 200000.5);
 295   fListTrig->Add(fHistCvlnTriggersPerRun);
 296
 297   fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 33000, 167000.5, 200000.5);
 298   fListTrig->Add(fHistMBTriggersPerRun);
 299
 300   fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 33000, 167000.5, 200000.5);
 301   fListTrig->Add(fHistCentralTriggersPerRun);
 302
 303   fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 33000, 167000.5, 200000.5);
 304   fListTrig->Add(fHistSemiCentralTriggersPerRun);
 305
 306   fListHist = new TList();
 307   fListHist ->SetOwner();
 308
 309   TString CutNameJPsi[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Two good tracks",
 310                                 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
 311   fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",13,0.5,13.5);
 312   for (Int_t i = 0; i<13; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
 313   fListHist->Add(fHistNeventsJPsi);
 314
 315   fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
 316   fListHist->Add(fHistTPCsignalJPsi);
 317
 318   fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
 319   fListHist->Add(fHistDiLeptonPtJPsi);
 320
 321   fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
 322   fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
 323   fListHist->Add(fHistDiElectronMass);
 324
 325   fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
 326   fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
 327   fListHist->Add(fHistDiMuonMass);
 328
 329   fHistDiLeptonMass = new TH1D("fHistDiLeptonMass","Invariant mass of J/#psi candidates",130,2.1,6.0);
 330   fHistDiLeptonMass->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}) (GeV/c)");
 331   fListHist->Add(fHistDiLeptonMass);
 332
 333   TString CutNamePsi2s[14] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Four good tracks",
 334                                 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
 335
 336   fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",14,0.5,14.5);
 337   for (Int_t i = 0; i<14; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
 338   fListHist->Add(fHistNeventsPsi2s);
 339
 340   fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
 341   fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
 342   fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 343   fListHist->Add(fHistPsi2sMassVsPt);
 344
 345   fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
 346   fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
 347   fListHist->Add(fHistPsi2sMassCoherent);
 348
 349   fListSystematics = new TList();
 350   fListSystematics->SetOwner();
 351   fListSystematics->SetName("fListSystematics");
 352   fListHist->Add(fListSystematics);
 353   InitSystematics();
 354
 355
 356   PostData(1, fJPsiTree);
 357   PostData(2, fPsi2sTree);
 358   PostData(3, fListTrig);
 359   PostData(4, fListHist);
 360
 361 }//UserCreateOutputObjects
 362
 363 //_____________________________________________________________________________
 364 void AliAnalysisTaskUpcPsi2s::InitSystematics()
 365 {
 366
 367 fListJPsiLoose = new TList();
 368 fListJPsiLoose->SetOwner();
 369 fListJPsiLoose->SetName("JPsiLoose");
 370 fListSystematics->Add(fListJPsiLoose);
 371
 372 TH1D *fHistJPsiNClusLoose = new TH1D("JPsiNClusLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 373 fListJPsiLoose->Add(fHistJPsiNClusLoose);
 374
 375 TH1D *fHistJPsiChi2Loose = new TH1D("JPsiChi2Loose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 376 fListJPsiLoose->Add(fHistJPsiChi2Loose);
 377
 378 TH1D *fHistJPsiDCAzLoose = new TH1D("JPsiDCAzLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 379 fListJPsiLoose->Add(fHistJPsiDCAzLoose);
 380
 381 TH1D *fHistJPsiDCAxyLoose = new TH1D("JPsiDCAxyLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 382 fListJPsiLoose->Add(fHistJPsiDCAxyLoose);
 383
 384 fListJPsiTight = new TList();
 385 fListJPsiTight->SetOwner();
 386 fListJPsiTight->SetName("JPsiTight");
 387 fListSystematics->Add(fListJPsiTight);
 388
 389 TH1D *fHistJPsiNClusTight = new TH1D("JPsiNClusTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
 390 fListJPsiTight->Add(fHistJPsiNClusTight);
 391
 392 TH1D *fHistJPsiChi2Tight = new TH1D("JPsiChi2Tight","Invariant mass of J/#psi candidates",130,2.1,6.0);
 393 fListJPsiTight->Add(fHistJPsiChi2Tight);
 394
 395 TH1D *fHistJPsiDCAzTight = new TH1D("JPsiDCAzTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
 396 fListJPsiTight->Add(fHistJPsiDCAzTight);
 397
 398 TH1D *fHistJPsiDCAxyTight = new TH1D("JPsiDCAxyTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
 399 fListJPsiTight->Add(fHistJPsiDCAxyTight);
 400
 401
 402 }
 403
 404 //_____________________________________________________________________________
 405 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
 406 {
 407
 408   //cout<<"#################### Next event ##################"<<endl;
 409
 410   if( fType == 0 ){
 411         RunESDtrig();
 412         if(fRunHist) RunESDhist();
 413         if(fRunTree) RunESDtree();
 414         }
 415
 416   if( fType == 1 ){
 417         RunAODtrig();
 418         if(fRunHist) RunAODhist();
 419         if(fRunTree) RunAODtree();
 420         }
 421
 422 }//UserExec
 423 //_____________________________________________________________________________
 424 void AliAnalysisTaskUpcPsi2s::RunAODtrig()
 425 {
 426
 427   //input event
 428   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 429   if(!aod) return;
 430
 431   fRunNum = aod ->GetRunNumber();
 432   //Trigger
 433   TString trigger = aod->GetFiredTriggerClasses();
 434
 435   if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
 436   if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
 437   if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
 438
 439   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
 440   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
 441
 442   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 443   if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
 444
 445   //MB, Central and SemiCentral triggers
 446   AliCentrality *centrality = aod->GetCentrality();
 447   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 448
 449   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
 450   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
 451
 452   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
 453
 454   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
 455
 456   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
 457
 458 PostData(3, fListTrig);
 459
 460 }
 461 //_____________________________________________________________________________
 462 void AliAnalysisTaskUpcPsi2s::RunAODhist()
 463 {
 464
 465   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 466
 467   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 468   Double_t muonMass = partMuon->Mass();
 469
 470   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 471   Double_t electronMass = partElectron->Mass();
 472
 473   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 474   Double_t pionMass = partPion->Mass();
 475
 476   //input event
 477   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 478   if(!aod) return;
 479
 480   cout<<"Event number: "<<((TTree*) GetInputData(0))->GetTree()->GetReadEntry()<<endl;
 481
 482   fHistNeventsJPsi->Fill(1);
 483   fHistNeventsPsi2s->Fill(1);
 484
 485   //Trigger
 486   TString trigger = aod->GetFiredTriggerClasses();
 487
 488   if(!isMC && !trigger.Contains("CCUP") ) return;
 489
 490   fHistNeventsJPsi->Fill(2);
 491   fHistNeventsPsi2s->Fill(2);
 492
 493   //primary vertex
 494   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 495   fVtxContrib = fAODVertex->GetNContributors();
 496   if(fVtxContrib < 2) return;
 497
 498   fHistNeventsJPsi->Fill(3);
 499   fHistNeventsPsi2s->Fill(3);
 500
 501   //VZERO, ZDC
 502   AliAODVZERO *fV0data = aod ->GetVZEROData();
 503   AliAODZDC *fZDCdata = aod->GetZDCData();
 504
 505   fV0Adecision = fV0data->GetV0ADecision();
 506   fV0Cdecision = fV0data->GetV0CDecision();
 507   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
 508
 509   fHistNeventsJPsi->Fill(4);
 510   fHistNeventsPsi2s->Fill(4);
 511
 512   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 513   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 514
 515   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
 516
 517   fHistNeventsJPsi->Fill(5);
 518   fHistNeventsPsi2s->Fill(5);
 519
 520   //Systematics - cut variation
 521   if(fRunSystematics) RunAODsystematics(aod);
 522
 523   //Two tracks loop
 524   Int_t nGoodTracks = 0;
 525   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 526
 527   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
 528   Short_t qLepton[4], qPion[4];
 529   UInt_t nLepton=0, nPion=0, nHighPt=0;
 530   Double_t fRecTPCsignal[5], fRecTPCsignalDist;
 531   Int_t fChannel = 0;
 532   Int_t mass[3]={-1,-1,-1};
 533
 534
 535   //Four track loop
 536   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 537     AliAODTrack *trk = aod->GetTrack(itr);
 538     if( !trk ) continue;
 539     if(!(trk->TestFilterBit(1<<0))) continue;
 540
 541       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 542       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 543       if(trk->GetTPCNcls() < 50)continue;
 544       if(trk->Chi2perNDF() > 4)continue;
 545       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 546       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 547       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 548       delete trk_clone;
 549
 550       if(TMath::Abs(dca[1]) > 2) continue;
 551
 552       TrackIndex[nGoodTracks] = itr;
 553       nGoodTracks++;
 554
 555       if(nGoodTracks > 4) break;
 556   }//Track loop
 557
 558   nLepton=0; nPion=0; nHighPt=0;
 559   mass[0]= -1; mass[1]= -1, mass[2]= -1;
 560
 561   if(nGoodTracks == 4){
 562           fHistNeventsPsi2s->Fill(6);
 563           for(Int_t i=0; i<4; i++){
 564                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 565
 566                 if(trk->Pt() > 1){
 567                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 568                         qLepton[nLepton] = trk->Charge();
 569                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 570                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 571                                         mass[nLepton] = 0;
 572                                         }
 573                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 574                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 575                                         mass[nLepton] = 1;
 576                                         }
 577                         nLepton++;
 578                         }
 579                 else{
 580                         qPion[nPion] = trk->Charge();
 581                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
 582                         nPion++;
 583                         }
 584
 585                 if(nLepton > 2 || nPion > 2) break;
 586                 }
 587         if((nLepton == 2) && (nPion == 2)){
 588                 fHistNeventsPsi2s->Fill(7);
 589                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(8);
 590                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(9);
 591                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
 592                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
 593                         fHistNeventsPsi2s->Fill(11);
 594                         if(mass[0] == mass[1]) {
 595                                 fHistNeventsPsi2s->Fill(12);
 596                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
 597                                 vDilepton = vLepton[0]+vLepton[1];
 598                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
 599                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
 600                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(13);
 601                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(14);
 602                                 }
 603                         }
 604                 }
 605   }
 606
 607   nGoodTracks = 0;
 608   //Two track loop
 609   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 610     AliAODTrack *trk = aod->GetTrack(itr);
 611     if( !trk ) continue;
 612     if(!(trk->TestFilterBit(1<<0))) continue;
 613
 614       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 615       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 616       if(trk->GetTPCNcls() < 70)continue;
 617       if(trk->Chi2perNDF() > 4)continue;
 618       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 619       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 620       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 621       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 622       delete trk_clone;
 623       if(TMath::Abs(dca[1]) > 2) continue;
 624       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 625       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
 626
 627       TrackIndex[nGoodTracks] = itr;
 628       nGoodTracks++;
 629
 630       if(nGoodTracks > 2) break;
 631   }//Track loop
 632
 633    nLepton=0; nPion=0; nHighPt=0;
 634   mass[0]= -1; mass[1]= -1, mass[2]= -1;
 635
 636   if(nGoodTracks == 2){
 637           fHistNeventsJPsi->Fill(6);
 638           for(Int_t i=0; i<2; i++){
 639                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 640                 if(trk->Pt() > 1) nHighPt++;
 641                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 642                 qLepton[nLepton] = trk->Charge();
 643                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 644                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 645                                 mass[nLepton] = 0;
 646                                 }
 647                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 648                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 649                                 mass[nLepton] = 1;
 650                                 }
 651                 nLepton++;
 652                 }
 653         if(nLepton == 2){
 654                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(7);
 655                 if(qLepton[0]*qLepton[1] < 0){
 656                         fHistNeventsJPsi->Fill(8);
 657                         if(nHighPt > 0){
 658                                 fHistNeventsJPsi->Fill(9);
 659                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
 660                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(10);
 661                                 if(mass[0] != -1 && mass[1] != -1) {
 662                                         fHistNeventsJPsi->Fill(11);
 663                                         vCandidate = vLepton[0]+vLepton[1];
 664                                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
 665                                         if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
 666                                         else {
 667                                                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
 668                                                 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
 669                                                 }
 670                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
 671                                         if(fChannel == -1) {
 672                                                 fHistDiMuonMass->Fill(vCandidate.M());
 673                                                 if(vCandidate.Pt()<0.15)fHistDiLeptonMass->Fill(vCandidate.M());
 674                                                 fHistNeventsJPsi->Fill(12);
 675                                                 }
 676                                         if(fChannel == 1) {
 677                                                 fHistDiElectronMass->Fill(vCandidate.M());
 678                                                 if(vCandidate.Pt()<0.3)fHistDiLeptonMass->Fill(vCandidate.M());
 679                                                 fHistNeventsJPsi->Fill(13);
 680                                                 }
 681                                         }
 682                                 }
 683                         }
 684                 }
 685   }
 686
 687
 688   PostData(4, fListHist);
 689
 690 }
 691
 692 //_____________________________________________________________________________
 693 void AliAnalysisTaskUpcPsi2s::RunAODsystematics(AliAODEvent* aod)
 694 {
 695
 696   Double_t fJPsiSels[4];
 697
 698   fJPsiSels[0] =   70; //min number of TPC clusters
 699   fJPsiSels[1] =   4; //chi2
 700   fJPsiSels[2] =   2; //DCAz
 701   fJPsiSels[3] =   1; // DCAxy 1x
 702
 703   Double_t fJPsiSelsMid[4];
 704
 705   fJPsiSelsMid[0] =   70; //min number of TPC clusters
 706   fJPsiSelsMid[1] =   4; //chi2
 707   fJPsiSelsMid[2] =   2; //DCAz
 708   fJPsiSelsMid[3] =   1; // DCAxy 1x
 709
 710   Double_t fJPsiSelsLoose[4];
 711
 712   fJPsiSelsLoose[0] =   60; //min number of TPC clusters
 713   fJPsiSelsLoose[1] =   5; //chi2
 714   fJPsiSelsLoose[2] =   3; //DCAz
 715   fJPsiSelsLoose[3] =   2; // DCAxy 2x
 716
 717   Double_t fJPsiSelsTight[4];
 718
 719   fJPsiSelsTight[0] =   80; //min number of TPC clusters
 720   fJPsiSelsTight[1] =   3.5; //chi2
 721   fJPsiSelsTight[2] =   1; //DCAz
 722   fJPsiSelsTight[3] =   0.5; // DCAxy 0.5x
 723
 724   Int_t nGoodTracks = 0;
 725   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 726
 727   TLorentzVector vLepton[4], vCandidate, vDilepton;
 728   Short_t qLepton[4];
 729   UInt_t nLepton=0, nHighPt=0;
 730   Double_t fRecTPCsignal[5], fRecTPCsignalDist;
 731   Int_t fChannel = 0;
 732
 733   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 734
 735   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 736
 737   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 738   Double_t muonMass = partMuon->Mass();
 739
 740   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 741   Double_t electronMass = partElectron->Mass();
 742
 743  // TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 744  // Double_t pionMass = partPion->Mass();
 745
 746
 747 for(Int_t i=0; i<4; i++){
 748           cout<<"Loose sytematics, cut"<<i<<endl;
 749           for(Int_t j=0; j<4; j++){
 750                   if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
 751                   else fJPsiSels[j] = fJPsiSelsMid[j];
 752           }
 753
 754   //Two track loop
 755   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 756     AliAODTrack *trk = aod->GetTrack(itr);
 757     if( !trk ) continue;
 758     if(!(trk->TestFilterBit(1<<0))) continue;
 759
 760       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 761       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 762       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) 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;
 766       delete trk_clone;
 767       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 768
 769       if(trk->GetTPCNcls() < fJPsiSels[0])continue;
 770       if(trk->Chi2perNDF() > fJPsiSels[1])continue;
 771       if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
 772       if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
 773
 774       TrackIndex[nGoodTracks] = itr;
 775       nGoodTracks++;
 776
 777       if(nGoodTracks > 2) break;
 778   }//Track loop
 779
 780   Int_t mass[3]={-1,-1,-1};
 781   fChannel = 0;
 782   nLepton=0; nHighPt=0;
 783
 784   if(nGoodTracks == 2){
 785           fHistNeventsJPsi->Fill(6);
 786           for(Int_t k=0; k<2; k++){
 787                 AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);
 788                 if(trk->Pt() > 1) nHighPt++;
 789                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 790                 qLepton[nLepton] = trk->Charge();
 791                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 792                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 793                                 mass[nLepton] = 0;
 794                                 }
 795                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 796                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 797                                 mass[nLepton] = 1;
 798                                 }
 799                 nLepton++;
 800                 }
 801         if(nLepton == 2){
 802                 if(qLepton[0]*qLepton[1] < 0 && nHighPt > 0 && (mass[0]!=-1 || mass[1]!=-1)){
 803                         vCandidate = vLepton[0]+vLepton[1];
 804                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
 805                         if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
 806                         else {
 807                                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
 808                                 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
 809                                 }
 810                         if(fChannel == -1 && vCandidate.Pt()<0.15) ((TH1D*)(fListJPsiLoose->At(i)))->Fill(vCandidate.M());
 811                         if(fChannel == 1 && vCandidate.Pt()<0.3) ((TH1D*)(fListJPsiLoose->At(i)))->Fill(vCandidate.M());
 812                         }
 813                 }
 814   }
 815 }//loose cuts
 816
 817 for(Int_t i=0; i<4; i++){
 818           cout<<"Tight sytematics, cut"<<i<<endl;
 819           for(Int_t j=0; j<4; j++){
 820                   if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
 821                   else fJPsiSels[j] = fJPsiSelsMid[j];
 822           }
 823
 824   //Two track loop
 825   nGoodTracks = 0;
 826   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 827     AliAODTrack *trk = aod->GetTrack(itr);
 828     if( !trk ) continue;
 829     if(!(trk->TestFilterBit(1<<0))) continue;
 830
 831       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 832       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 833       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 834       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 835       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 836       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 837       delete trk_clone;
 838       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 839
 840       if(trk->GetTPCNcls() < fJPsiSels[0])continue;
 841       if(trk->Chi2perNDF() > fJPsiSels[1])continue;
 842       if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
 843       if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
 844
 845       TrackIndex[nGoodTracks] = itr;
 846       nGoodTracks++;
 847
 848       if(nGoodTracks > 2) break;
 849   }//Track loop
 850
 851   Int_t mass[3]={-1,-1,-1};
 852   fChannel = 0;
 853   nLepton=0; nHighPt=0;
 854
 855   if(nGoodTracks == 2){
 856           fHistNeventsJPsi->Fill(6);
 857           for(Int_t k=0; k<2; k++){
 858                 AliAODTrack *trk = aod->GetTrack(TrackIndex[k]);
 859                 if(trk->Pt() > 1) nHighPt++;
 860                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 861                 qLepton[nLepton] = trk->Charge();
 862                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 863                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 864                                 mass[nLepton] = 0;
 865                                 }
 866                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 867                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 868                                 mass[nLepton] = 1;
 869                                 }
 870                 nLepton++;
 871                 }
 872         if(nLepton == 2){
 873                 if(qLepton[0]*qLepton[1] < 0 && nHighPt > 0 && (mass[0]!=-1 || mass[1]!=-1)){
 874                         vCandidate = vLepton[0]+vLepton[1];
 875                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
 876                         if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
 877                         else {
 878                                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
 879                                 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
 880                                 }
 881                         if(fChannel == -1 && vCandidate.Pt()<0.15) ((TH1D*)(fListJPsiTight->At(i)))->Fill(vCandidate.M());
 882                         if(fChannel == 1 && vCandidate.Pt()<0.3) ((TH1D*)(fListJPsiTight->At(i)))->Fill(vCandidate.M());
 883                         }
 884                 }
 885   }
 886 }//tight cuts
 887
 888
 889 }
 890 //_____________________________________________________________________________
 891 void AliAnalysisTaskUpcPsi2s::RunAODtree()
 892 {
 893   //input event
 894   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 895   if(!aod) return;
 896
 897   if(isMC) RunAODMC(aod);
 898
 899   //input data
 900   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
 901   fDataFilnam->Clear();
 902   fDataFilnam->SetString(filnam);
 903   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
 904   fRunNum = aod ->GetRunNumber();
 905
 906   //Trigger
 907   TString trigger = aod->GetFiredTriggerClasses();
 908
 909   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
 910   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
 911   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
 912
 913   Bool_t isTriggered = kFALSE;
 914   for(Int_t i=0; i<ntrg; i++) {
 915     if( fTrigger[i] ) isTriggered = kTRUE;
 916   }
 917   if(!isMC && !isTriggered ) return;
 918
 919   //trigger inputs
 920   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
 921   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 922
 923   //Event identification
 924   fPerNum = aod ->GetPeriodNumber();
 925   fOrbNum = aod ->GetOrbitNumber();
 926   fBCrossNum = aod ->GetBunchCrossNumber();
 927
 928   //primary vertex
 929   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 930   fVtxContrib = fAODVertex->GetNContributors();
 931   fVtxPos[0] = fAODVertex->GetX();
 932   fVtxPos[1] = fAODVertex->GetY();
 933   fVtxPos[2] = fAODVertex->GetZ();
 934   Double_t CovMatx[6];
 935   fAODVertex->GetCovarianceMatrix(CovMatx);
 936   fVtxErr[0] = CovMatx[0];
 937   fVtxErr[1] = CovMatx[1];
 938   fVtxErr[2] = CovMatx[2];
 939   fVtxChi2 = fAODVertex->GetChi2();
 940   fVtxNDF = fAODVertex->GetNDF();
 941
 942   //Tracklets
 943   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
 944
 945   //VZERO, ZDC
 946   AliAODVZERO *fV0data = aod ->GetVZEROData();
 947   AliAODZDC *fZDCdata = aod->GetZDCData();
 948
 949   fV0Adecision = fV0data->GetV0ADecision();
 950   fV0Cdecision = fV0data->GetV0CDecision();
 951   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 952   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 953
 954   fNLooseTracks = 0;
 955
 956   //Track loop - loose cuts
 957   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 958     AliAODTrack *trk = aod->GetTrack(itr);
 959     if( !trk ) continue;
 960     if(!(trk->TestFilterBit(1<<0))) continue;
 961
 962       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 963       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 964       if(trk->GetTPCNcls() < 20)continue;
 965       fNLooseTracks++;
 966   }//Track loop -loose cuts
 967
 968   Int_t nGoodTracks=0;
 969   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 970
 971   //Two track loop
 972   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 973     AliAODTrack *trk = aod->GetTrack(itr);
 974     if( !trk ) continue;
 975     if(!(trk->TestFilterBit(1<<0))) continue;
 976
 977       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 978       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 979       if(trk->GetTPCNcls() < 70)continue;
 980       if(trk->Chi2perNDF() > 4)continue;
 981       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 982       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 983       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 984       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 985       delete trk_clone;
 986       if(TMath::Abs(dca[1]) > 2) continue;
 987       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 988       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
 989
 990
 991       TrackIndex[nGoodTracks] = itr;
 992       nGoodTracks++;
 993
 994       if(nGoodTracks > 2) break;
 995   }//Track loop
 996
 997   fJPsiAODTracks->Clear("C");
 998   if(nGoodTracks == 2){
 999
1000           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1001           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1002           Double_t muonMass = partMuon->Mass();
1003           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1004           Double_t electronMass = partElectron->Mass();
1005           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1006           Double_t pionMass = partPion->Mass();
1007
1008           Double_t KFcov[21];
1009           Double_t KFpar[6];
1010           Double_t KFmass = pionMass;
1011           Double_t fRecTPCsignal;
1012           AliKFParticle *KFpart[2];
1013           AliKFVertex *KFvtx = new AliKFVertex();
1014           KFvtx->SetField(aod->GetMagneticField());
1015
1016           for(Int_t i=0; i<2; i++){
1017                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
1018
1019                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1020                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1021                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1022                 delete trk_clone;
1023
1024                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
1025                 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
1026
1027                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1028                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1029                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1030
1031                 trk->GetPosition(KFpar);
1032                 trk->PxPyPz(KFpar+3);
1033                 trk->GetCovMatrix(KFcov);
1034
1035                 if(trk->Pt() > 1){
1036                         fRecTPCsignal = trk->GetTPCsignal();
1037                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
1038                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
1039                         }
1040                 else KFmass = pionMass;
1041
1042                 KFpart[i] = new AliKFParticle();
1043                 KFpart[i]->SetField(aod->GetMagneticField());
1044                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
1045                 KFvtx->AddDaughter(*KFpart[i]);
1046
1047
1048                 Double_t pos[3]={0,0,0};
1049                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
1050                 parTrk->CopyFromVTrack((AliVTrack*) trk);
1051                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
1052                 else {
1053                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1054                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1055                      }
1056                 delete parTrk;
1057                 }
1058   fKfVtxPos[0]= KFvtx->GetX();
1059   fKfVtxPos[1]= KFvtx->GetY();
1060   fKfVtxPos[2]= KFvtx->GetZ();
1061   for(UInt_t i=0; i<2; i++)delete KFpart[i];
1062   delete KFvtx;
1063
1064   if(!isMC) fJPsiTree ->Fill();
1065   }
1066
1067    nGoodTracks = 0;
1068    //Four track loop
1069   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1070     AliAODTrack *trk = aod->GetTrack(itr);
1071     if( !trk ) continue;
1072     if(!(trk->TestFilterBit(1<<0))) continue;
1073
1074       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
1075       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
1076       if(trk->GetTPCNcls() < 50)continue;
1077       if(trk->Chi2perNDF() > 4)continue;
1078       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1079       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1080       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1081       delete trk_clone;
1082       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1083       if(TMath::Abs(dca[1]) > 2) continue;
1084       Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1085       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
1086
1087       TrackIndex[nGoodTracks] = itr;
1088       nGoodTracks++;
1089
1090       if(nGoodTracks > 4) break;
1091   }//Track loop
1092
1093   fPsi2sAODTracks->Clear("C");
1094   if(nGoodTracks == 4){
1095
1096           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1097           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1098           Double_t muonMass = partMuon->Mass();
1099           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1100           Double_t electronMass = partElectron->Mass();
1101           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1102           Double_t pionMass = partPion->Mass();
1103
1104           Double_t KFcov[21];
1105           Double_t KFpar[6];
1106           Double_t KFmass = pionMass;
1107           Double_t fRecTPCsignal;
1108           AliKFParticle *KFpart[4];
1109           AliKFVertex *KFvtx = new AliKFVertex();
1110           KFvtx->SetField(aod->GetMagneticField());
1111
1112           for(Int_t i=0; i<4; i++){
1113                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
1114
1115                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1116                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1117                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1118                 delete trk_clone;
1119
1120                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
1121                 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
1122
1123
1124                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1125                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1126                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1127
1128                 trk->GetPosition(KFpar);
1129                 trk->PxPyPz(KFpar+3);
1130                 trk->GetCovMatrix(KFcov);
1131
1132                 if(trk->Pt() > 1){
1133                         fRecTPCsignal = trk->GetTPCsignal();
1134                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
1135                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
1136                         }
1137                 else KFmass = pionMass;
1138
1139                 KFpart[i] = new AliKFParticle();
1140                 KFpart[i]->SetField(aod->GetMagneticField());
1141                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
1142                 KFvtx->AddDaughter(*KFpart[i]);
1143
1144                 Double_t pos[3]={0,0,0};
1145                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
1146                 parTrk->CopyFromVTrack((AliVTrack*) trk);
1147                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
1148                 else {
1149                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1150                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1151                      }
1152                 delete parTrk;
1153                 }
1154   fKfVtxPos[0]= KFvtx->GetX();
1155   fKfVtxPos[1]= KFvtx->GetY();
1156   fKfVtxPos[2]= KFvtx->GetZ();
1157   for(UInt_t i=0; i<4; i++)delete KFpart[i];
1158   delete KFvtx;
1159   if(!isMC) fPsi2sTree ->Fill();
1160   }
1161
1162   if(isMC){
1163         fJPsiTree ->Fill();
1164         fPsi2sTree ->Fill();
1165   }
1166
1167   PostData(1, fJPsiTree);
1168   PostData(2, fPsi2sTree);
1169
1170 }//RunAOD
1171
1172
1173 //_____________________________________________________________________________
1174 void AliAnalysisTaskUpcPsi2s::RunAODMC(AliAODEvent *aod)
1175 {
1176
1177   fGenPart->Clear("C");
1178
1179   TClonesArray *arrayMC = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1180   if(!arrayMC) return;
1181
1182   Int_t nmc=0;
1183   //loop over mc particles
1184   for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
1185     AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(imc);
1186     if(!mcPart) continue;
1187
1188     if(mcPart->GetMother() >= 0) continue;
1189
1190     TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1191     part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1192     part->SetPdgCode(mcPart->GetPdgCode());
1193     part->SetUniqueID(imc);
1194   }//loop over mc particles
1195
1196 }//RunAODMC
1197
1198
1199 //_____________________________________________________________________________
1200 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
1201 {
1202
1203   //input event
1204   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1205   if(!esd) return;
1206
1207   fRunNum = esd ->GetRunNumber();
1208   //Trigger
1209   TString trigger = esd->GetFiredTriggerClasses();
1210
1211   if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
1212   if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
1213   if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
1214
1215   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
1216   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
1217
1218   if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
1219
1220    //MB, Central and SemiCentral triggers
1221   AliCentrality *centrality = esd->GetCentrality();
1222   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1223
1224   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
1225   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
1226
1227   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
1228
1229   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
1230
1231   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
1232
1233
1234 PostData(3, fListTrig);
1235
1236 }
1237 //_____________________________________________________________________________
1238 void AliAnalysisTaskUpcPsi2s::RunESDhist()
1239 {
1240
1241   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1242
1243   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1244   Double_t muonMass = partMuon->Mass();
1245
1246   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1247   Double_t electronMass = partElectron->Mass();
1248
1249   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1250   Double_t pionMass = partPion->Mass();
1251
1252   //input event
1253   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1254   if(!esd) return;
1255
1256   fHistNeventsJPsi->Fill(1);
1257   fHistNeventsPsi2s->Fill(1);
1258
1259   //Trigger
1260   TString trigger = esd->GetFiredTriggerClasses();
1261
1262   if(!isMC && !trigger.Contains("CCUP") ) return;
1263
1264   fHistNeventsJPsi->Fill(2);
1265   fHistNeventsPsi2s->Fill(2);
1266
1267   //primary vertex
1268   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1269   fVtxContrib = fESDVertex->GetNContributors();
1270   if(fVtxContrib < 2) return;
1271
1272   fHistNeventsJPsi->Fill(3);
1273   fHistNeventsPsi2s->Fill(3);
1274
1275   //VZERO, ZDC
1276   AliESDVZERO *fV0data = esd->GetVZEROData();
1277   AliESDZDC *fZDCdata = esd->GetESDZDC();
1278
1279   fV0Adecision = fV0data->GetV0ADecision();
1280   fV0Cdecision = fV0data->GetV0CDecision();
1281   if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
1282
1283   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1284   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1285   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
1286
1287   fHistNeventsJPsi->Fill(4);
1288   fHistNeventsPsi2s->Fill(4);
1289
1290    Int_t nGoodTracks=0;
1291   //Two tracks loop
1292   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1293
1294   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1295   Short_t qLepton[4], qPion[4];
1296   UInt_t nLepton=0, nPion=0, nHighPt=0;
1297   Double_t fRecTPCsignal[5];
1298   Int_t mass[3]={-1,-1,-1};
1299
1300  //Two Track loop
1301   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1302     AliESDtrack *trk = esd->GetTrack(itr);
1303     if( !trk ) continue;
1304
1305       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1306       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1307       if(trk->GetTPCNcls() < 70)continue;
1308       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1309       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1310       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1311       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1312       trk->GetImpactParameters(dca[0],dca[1]);
1313       if(TMath::Abs(dca[1]) > 2) continue;
1314       if(TMath::Abs(dca[1]) > 0.2) continue;
1315
1316       TrackIndex[nGoodTracks] = itr;
1317       nGoodTracks++;
1318       if(nGoodTracks > 2) break;
1319   }//Track loop
1320
1321
1322   if(nGoodTracks == 2){
1323           fHistNeventsJPsi->Fill(5);
1324           for(Int_t i=0; i<2; i++){
1325                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1326                 if(trk->Pt() > 1) nHighPt++;
1327                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1328                 qLepton[nLepton] = trk->Charge();
1329                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1330                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1331                                 mass[nLepton] = 0;
1332                                 }
1333                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1334                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1335                                 mass[nLepton] = 1;
1336                                 }
1337                 nLepton++;
1338                 }
1339         if(nLepton == 2){
1340                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
1341                 if(qLepton[0]*qLepton[1] < 0){
1342                         fHistNeventsJPsi->Fill(7);
1343                         if(nHighPt > 0){
1344                                 fHistNeventsJPsi->Fill(8);
1345                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
1346                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
1347                                 if(mass[0] == mass[1] && mass[0] != -1) {
1348                                         fHistNeventsJPsi->Fill(10);
1349                                         vCandidate = vLepton[0]+vLepton[1];
1350                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
1351                                         if(mass[0] == 0) {
1352                                                 fHistDiMuonMass->Fill(vCandidate.M());
1353                                                 fHistNeventsJPsi->Fill(11);
1354                                                 }
1355                                         if(mass[0] == 1) {
1356                                                 fHistDiElectronMass->Fill(vCandidate.M());
1357                                                 fHistNeventsJPsi->Fill(12);
1358                                                 }
1359                                         }
1360                                 }
1361                         }
1362                 }
1363   }
1364   nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1365   mass[0]= -1; mass[1]= -1, mass[2]= -1;
1366
1367     //Four Track loop
1368   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1369     AliESDtrack *trk = esd->GetTrack(itr);
1370     if( !trk ) continue;
1371
1372       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1373       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1374       if(trk->GetTPCNcls() < 50)continue;
1375       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1376       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1377       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1378       trk->GetImpactParameters(dca[0],dca[1]);
1379       if(TMath::Abs(dca[1]) > 2) continue;
1380
1381       TrackIndex[nGoodTracks] = itr;
1382       nGoodTracks++;
1383       if(nGoodTracks > 4) break;
1384   }//Track loop
1385
1386   if(nGoodTracks == 4){
1387           fHistNeventsPsi2s->Fill(5);
1388           for(Int_t i=0; i<4; i++){
1389                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1390
1391                 if(trk->Pt() > 1){
1392                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1393                         qLepton[nLepton] = trk->Charge();
1394                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1395                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1396                                         mass[nLepton] = 0;
1397                                         }
1398                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1399                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1400                                         mass[nLepton] = 1;
1401                                         }
1402                         nLepton++;
1403                         }
1404                 else{
1405                         qPion[nPion] = trk->Charge();
1406                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1407                         nPion++;
1408                         }
1409
1410                 if(nLepton > 2 || nPion > 2) break;
1411                 }
1412         if((nLepton == 2) && (nPion == 2)){
1413                 fHistNeventsPsi2s->Fill(6);
1414                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1415                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1416                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1417                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1418                         fHistNeventsPsi2s->Fill(10);
1419                         if(mass[0] == mass[1]) {
1420                                 fHistNeventsPsi2s->Fill(11);
1421                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1422                                 vDilepton = vLepton[0]+vLepton[1];
1423                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1424                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1425                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1426                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1427                                 }
1428                         }
1429                 }
1430   }
1431
1432   PostData(4, fListHist);
1433
1434 }
1435
1436 //_____________________________________________________________________________
1437 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1438 {
1439
1440   //input event
1441   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1442   if(!esd) return;
1443
1444   if(isMC) RunESDMC(esd);
1445
1446   //input data
1447   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1448   fDataFilnam->Clear();
1449   fDataFilnam->SetString(filnam);
1450   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1451   fRunNum = esd->GetRunNumber();
1452
1453    //Trigger
1454   TString trigger = esd->GetFiredTriggerClasses();
1455
1456   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1457   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
1458   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1459
1460   Bool_t isTriggered = kFALSE;
1461   for(Int_t i=0; i<ntrg; i++) {
1462     if( fTrigger[i] ) isTriggered = kTRUE;
1463   }
1464   if(!isMC && !isTriggered ) return;
1465
1466   //trigger inputs
1467   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1468   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1469
1470   //TOF trigger info (0OMU)
1471   fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1472   fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1473
1474   //Event identification
1475   fPerNum = esd->GetPeriodNumber();
1476   fOrbNum = esd->GetOrbitNumber();
1477   fBCrossNum = esd->GetBunchCrossNumber();
1478
1479   //primary vertex
1480   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1481   fVtxContrib = fESDVertex->GetNContributors();
1482   fVtxPos[0] = fESDVertex->GetX();
1483   fVtxPos[1] = fESDVertex->GetY();
1484   fVtxPos[2] = fESDVertex->GetZ();
1485   Double_t CovMatx[6];
1486   fESDVertex->GetCovarianceMatrix(CovMatx);
1487   fVtxErr[0] = CovMatx[0];
1488   fVtxErr[1] = CovMatx[1];
1489   fVtxErr[2] = CovMatx[2];
1490   fVtxChi2 = fESDVertex->GetChi2();
1491   fVtxNDF = fESDVertex->GetNDF();
1492
1493   //Tracklets
1494   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1495
1496   //VZERO, ZDC
1497   AliESDVZERO *fV0data = esd->GetVZEROData();
1498   AliESDZDC *fZDCdata = esd->GetESDZDC();
1499
1500   fV0Adecision = fV0data->GetV0ADecision();
1501   fV0Cdecision = fV0data->GetV0CDecision();
1502   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1503   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1504
1505   fNLooseTracks = 0;
1506
1507   //Track loop - loose cuts
1508   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1509     AliESDtrack *trk = esd->GetTrack(itr);
1510     if( !trk ) continue;
1511
1512       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1513       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1514       if(trk->GetTPCNcls() < 20)continue;
1515       fNLooseTracks++;
1516   }//Track loop -loose cuts
1517
1518   Int_t nGoodTracks=0;
1519   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1520
1521   //Two Track loop
1522   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1523     AliESDtrack *trk = esd->GetTrack(itr);
1524     if( !trk ) continue;
1525
1526       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1527       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1528       if(trk->GetTPCNcls() < 70)continue;
1529       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1530       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1531       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1532       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1533       trk->GetImpactParameters(dca[0],dca[1]);
1534       if(TMath::Abs(dca[1]) > 2) continue;
1535       if(TMath::Abs(dca[1]) > 0.2) continue;
1536
1537       TrackIndex[nGoodTracks] = itr;
1538       nGoodTracks++;
1539       if(nGoodTracks > 2) break;
1540   }//Track loop
1541
1542   fJPsiESDTracks->Clear("C");
1543   if(nGoodTracks == 2){
1544           for(Int_t i=0; i<2; i++){
1545                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1546
1547                 AliExternalTrackParam cParam;
1548                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1549
1550                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
1551
1552                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1553                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1554                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1555
1556                 Double_t pos[3]={0,0,0};
1557                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1558                 else {
1559                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1560                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1561                      }
1562                 }
1563   if(!isMC) fJPsiTree ->Fill();
1564   }
1565
1566   nGoodTracks = 0;
1567   //Four track loop
1568   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1569     AliESDtrack *trk = esd->GetTrack(itr);
1570     if( !trk ) continue;
1571
1572       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1573       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1574       if(trk->GetTPCNcls() < 50)continue;
1575       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1576       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1577       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1578       trk->GetImpactParameters(dca[0],dca[1]);
1579       if(TMath::Abs(dca[1]) > 2) continue;
1580
1581       TrackIndex[nGoodTracks] = itr;
1582       nGoodTracks++;
1583       if(nGoodTracks > 4) break;
1584   }//Track loop
1585
1586   fPsi2sESDTracks->Clear("C");
1587   if(nGoodTracks == 4){
1588           for(Int_t i=0; i<4; i++){
1589                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1590
1591                 AliExternalTrackParam cParam;
1592                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1593
1594                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1595
1596                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1597                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1598                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1599
1600                 Double_t pos[3]={0,0,0};
1601                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1602                 else {
1603                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1604                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1605                      }
1606                 }
1607   if(!isMC) fPsi2sTree ->Fill();
1608   }
1609
1610   if(isMC){
1611         fJPsiTree ->Fill();
1612         fPsi2sTree ->Fill();
1613   }
1614
1615   PostData(1, fJPsiTree);
1616   PostData(2, fPsi2sTree);
1617
1618 }//RunESD
1619
1620
1621 //_____________________________________________________________________________
1622 void AliAnalysisTaskUpcPsi2s::RunESDMC(AliESDEvent* esd)
1623 {
1624
1625   AliTriggerAnalysis *fTrigAna = new AliTriggerAnalysis();
1626   fTrigAna->SetAnalyzeMC(isMC);
1627
1628   if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) > 1) fTriggerInputsMC[0] = kTRUE;
1629   if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) < 2) fTriggerInputsMC[0] = kFALSE;
1630
1631   fTriggerInputsMC[1] = esd->GetHeader()->IsTriggerInputFired("0OMU");
1632   fTriggerInputsMC[2] = esd->GetHeader()->IsTriggerInputFired("0VBA");
1633   fTriggerInputsMC[3] = esd->GetHeader()->IsTriggerInputFired("0VBC");
1634
1635   fGenPart->Clear("C");
1636
1637   AliMCEvent *mc = MCEvent();
1638   if(!mc) return;
1639
1640   Int_t nmc = 0;
1641   //loop over mc particles
1642   for(Int_t imc=0; imc<mc->GetNumberOfTracks(); imc++) {
1643     AliMCParticle *mcPart = (AliMCParticle*) mc->GetTrack(imc);
1644     if(!mcPart) continue;
1645
1646     if(mcPart->GetMother() >= 0) continue;
1647
1648     TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1649     part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1650     part->SetPdgCode(mcPart->PdgCode());
1651     part->SetUniqueID(imc);
1652   }//loop over mc particles
1653
1654 }//RunESDMC
1655
1656
1657
1658 //_____________________________________________________________________________
1659 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
1660 {
1661
1662   cout<<"Analysis complete."<<endl;
1663 }//Terminate
1664