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),fSpdVtxContrib(0),
  72     fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
  73     fZDCAenergy(0),fZDCCenergy(0),fZDCAtime(0),fZDCCtime(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),fHistCint1TriggersPerRun(0),fHistC0tvxAndCint1TriggersPerRun(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),fListPsi2sLoose(0),fListPsi2sTight(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),fSpdVtxContrib(0),
  95     fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
  96     fZDCAenergy(0),fZDCCenergy(0),fZDCAtime(0),fZDCCtime(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),fHistCint1TriggersPerRun(0),fHistC0tvxAndCint1TriggersPerRun(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),fListPsi2sLoose(0),fListPsi2sTight(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         fPIDTPCMuon[i] = -666;
 128         fPIDTPCElectron[i] = -666;
 129         fPIDTPCPion[i] = -666;
 130         fPIDTPCKaon[i] = -666;
 131         fPIDTPCProton[i] = -666;
 132
 133         fPIDTOFMuon[i] = -666;
 134         fPIDTOFElectron[i] = -666;
 135         fPIDTOFPion[i] = -666;
 136         fPIDTOFKaon[i] = -666;
 137         fPIDTOFProton[i] = -666;
 138
 139         fTriggerInputsMC[i] = kFALSE;
 140
 141         fIsVtxContributor[i] = kFALSE;
 142         }
 143   for(Int_t i=0; i<3; i++){
 144         fVtxPos[i] = -666;
 145         fVtxErr[i] = -666;
 146         fKfVtxPos[i] = -666;
 147         fSpdVtxPos[i] = -666;
 148         }
 149
 150 }//Init
 151
 152 //_____________________________________________________________________________
 153 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
 154 {
 155   // Destructor
 156   if(fJPsiTree){
 157      delete fJPsiTree;
 158      fJPsiTree = 0x0;
 159   }
 160   if(fPsi2sTree){
 161      delete fPsi2sTree;
 162      fPsi2sTree = 0x0;
 163   }
 164   if(fListTrig){
 165      delete fListTrig;
 166      fListTrig = 0x0;
 167   }
 168   if(fListHist){
 169      delete fListHist;
 170      fListHist = 0x0;
 171   }
 172
 173 }//~AliAnalysisTaskUpcPsi2s
 174
 175
 176 //_____________________________________________________________________________
 177 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
 178 {
 179   //PID response
 180   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 181   AliInputEventHandler *inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 182   fPIDResponse = inputHandler->GetPIDResponse();
 183
 184   //input file
 185   fDataFilnam = new TObjString();
 186   fDataFilnam->SetString("");
 187
 188     //tracks
 189   fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
 190   fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
 191   fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
 192   fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
 193   fGenPart = new TClonesArray("TParticle", 1000);
 194
 195   //output tree with JPsi candidate events
 196   fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
 197   fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 198   fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 199   fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 200
 201   fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 202   fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 203   fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 204   fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 205   fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 206   fJPsiTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
 207   fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 208   fJPsiTree ->Branch("fSpdVtxContrib", &fSpdVtxContrib, "fSpdVtxContrib/I");
 209
 210   fJPsiTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
 211   fJPsiTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
 212   fJPsiTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[2]/D");
 213
 214   fJPsiTree ->Branch("fPIDTPCMuon", &fPIDTPCMuon[0], "fPIDTPCMuon[2]/D");
 215   fJPsiTree ->Branch("fPIDTPCElectron", &fPIDTPCElectron[0], "fPIDTPCElectron[2]/D");
 216   fJPsiTree ->Branch("fPIDTPCPion", &fPIDTPCPion[0], "fPIDTPCPion[2]/D");
 217   fJPsiTree ->Branch("fPIDTPCKaon", &fPIDTPCKaon[0], "fPIDTPCKaon[2]/D");
 218   fJPsiTree ->Branch("fPIDTPCProton", &fPIDTPCProton[0], "fPIDTPCProton[2]/D");
 219
 220   fJPsiTree ->Branch("fPIDTOFMuon", &fPIDTOFMuon[0], "fPIDTOFMuon[2]/D");
 221   fJPsiTree ->Branch("fPIDTOFElectron", &fPIDTOFElectron[0], "fPIDTOFElectron[2]/D");
 222   fJPsiTree ->Branch("fPIDTOFPion", &fPIDTOFPion[0], "fPIDTOFPion[2]/D");
 223   fJPsiTree ->Branch("fPIDTOFKaon", &fPIDTOFKaon[0], "fPIDTOFKaon[2]/D");
 224   fJPsiTree ->Branch("fPIDTOFProton", &fPIDTOFProton[0], "fPIDTOFProton[2]/D");
 225
 226   fJPsiTree ->Branch("fIsVtxContributor", &fIsVtxContributor[0], "fIsVtxContributor[2]/O");
 227
 228   fJPsiTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
 229   fJPsiTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
 230   fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
 231   fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
 232
 233   fJPsiTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
 234   fJPsiTree ->Branch("fSpdVtxPos", &fSpdVtxPos[0], "fSpdVtxPos[3]/D");
 235
 236   fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 237   fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 238   fJPsiTree ->Branch("fZDCAtime", &fZDCAtime, "fZDCAtime/D");
 239   fJPsiTree ->Branch("fZDCCtime", &fZDCCtime, "fZDCCtime/D");
 240   fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 241   fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 242   fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
 243   fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 244   fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 245   if( fType == 0 ) {
 246     fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
 247   }
 248   if( fType == 1 ) {
 249     fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
 250   }
 251   if(isMC) {
 252     fJPsiTree ->Branch("fGenPart", &fGenPart);
 253     fJPsiTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
 254   }
 255
 256
 257  //output tree with Psi2s candidate events
 258   fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
 259   fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 260   fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 261   fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 262
 263   fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 264   fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 265   fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 266   fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 267   fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 268   fPsi2sTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
 269   fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 270   fPsi2sTree ->Branch("fSpdVtxContrib", &fSpdVtxContrib, "fSpdVtxContrib/I");
 271
 272   fPsi2sTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
 273   fPsi2sTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
 274   fPsi2sTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
 275
 276   fPsi2sTree ->Branch("fPIDTPCMuon", &fPIDTPCMuon[0], "fPIDTPCMuon[4]/D");
 277   fPsi2sTree ->Branch("fPIDTPCElectron", &fPIDTPCElectron[0], "fPIDTPCElectron[4]/D");
 278   fPsi2sTree ->Branch("fPIDTPCPion", &fPIDTPCPion[0], "fPIDTPCPion[4]/D");
 279   fPsi2sTree ->Branch("fPIDTPCKaon", &fPIDTPCKaon[0], "fPIDTPCKaon[4]/D");
 280   fPsi2sTree ->Branch("fPIDTPCProton", &fPIDTPCProton[0], "fPIDTPCProton[4]/D");
 281
 282   fPsi2sTree ->Branch("fPIDTOFMuon", &fPIDTOFMuon[0], "fPIDTOFMuon[4]/D");
 283   fPsi2sTree ->Branch("fPIDTOFElectron", &fPIDTOFElectron[0], "fPIDTOFElectron[4]/D");
 284   fPsi2sTree ->Branch("fPIDTOFPion", &fPIDTOFPion[0], "fPIDTOFPion[4]/D");
 285   fPsi2sTree ->Branch("fPIDTOFKaon", &fPIDTOFKaon[0], "fPIDTOFKaon[4]/D");
 286   fPsi2sTree ->Branch("fPIDTOFProton", &fPIDTOFProton[0], "fPIDTOFProton[4]/D");
 287
 288   fPsi2sTree ->Branch("fIsVtxContributor", &fIsVtxContributor[0], "fIsVtxContributor[4]/O");
 289
 290   fPsi2sTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
 291   fPsi2sTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
 292   fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
 293   fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
 294
 295   fPsi2sTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
 296   fPsi2sTree ->Branch("fSpdVtxPos", &fSpdVtxPos[0], "fSpdVtxPos[3]/D");
 297
 298   fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 299   fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 300   fPsi2sTree ->Branch("fZDCAtime", &fZDCAtime, "fZDCAtime/D");
 301   fPsi2sTree ->Branch("fZDCCtime", &fZDCCtime, "fZDCCtime/D");
 302   fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 303   fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 304   fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
 305   fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 306   fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 307   if( fType == 0 ) {
 308     fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
 309   }
 310   if( fType == 1 ) {
 311     fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
 312   }
 313   if(isMC) {
 314     fPsi2sTree ->Branch("fGenPart", &fGenPart);
 315     fPsi2sTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
 316   }
 317
 318
 319   fListTrig = new TList();
 320   fListTrig ->SetOwner();
 321
 322   fHistCcup4TriggersPerRun = new TH1D("fHistCcup4TriggersPerRun", "fHistCcup4TriggersPerRun", 33000, 167000.5, 200000.5);
 323   fListTrig->Add(fHistCcup4TriggersPerRun);
 324
 325   fHistCcup7TriggersPerRun = new TH1D("fHistCcup7TriggersPerRun", "fHistCcup7TriggersPerRun", 33000, 167000.5, 200000.5);
 326   fListTrig->Add(fHistCcup7TriggersPerRun);
 327
 328   fHistCcup2TriggersPerRun = new TH1D("fHistCcup2TriggersPerRun", "fHistCcup2TriggersPerRun", 33000, 167000.5, 200000.5);
 329   fListTrig->Add(fHistCcup2TriggersPerRun);
 330
 331   fHistCint1TriggersPerRun = new TH1D("fHistCint1TriggersPerRun", "fHistCint1TriggersPerRun", 33000, 167000.5, 200000.5);
 332   fListTrig->Add(fHistCint1TriggersPerRun);
 333
 334   fHistC0tvxAndCint1TriggersPerRun = new TH1D("fHistC0tvxAndCint1TriggersPerRun", "fHistC0tvxAndCint1TriggersPerRun", 33000, 167000.5, 200000.5);
 335   fListTrig->Add(fHistC0tvxAndCint1TriggersPerRun);
 336
 337   fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 33000, 167000.5, 200000.5);
 338   fListTrig->Add(fHistZedTriggersPerRun);
 339
 340   fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 33000, 167000.5, 200000.5);
 341   fListTrig->Add(fHistCvlnTriggersPerRun);
 342
 343   fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 33000, 167000.5, 200000.5);
 344   fListTrig->Add(fHistMBTriggersPerRun);
 345
 346   fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 33000, 167000.5, 200000.5);
 347   fListTrig->Add(fHistCentralTriggersPerRun);
 348
 349   fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 33000, 167000.5, 200000.5);
 350   fListTrig->Add(fHistSemiCentralTriggersPerRun);
 351
 352   fListHist = new TList();
 353   fListHist ->SetOwner();
 354
 355   TString CutNameJPsi[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Two good tracks",
 356                                 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
 357   fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",13,0.5,13.5);
 358   for (Int_t i = 0; i<13; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
 359   fListHist->Add(fHistNeventsJPsi);
 360
 361   fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
 362   fListHist->Add(fHistTPCsignalJPsi);
 363
 364   fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
 365   fListHist->Add(fHistDiLeptonPtJPsi);
 366
 367   fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
 368   fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
 369   fListHist->Add(fHistDiElectronMass);
 370
 371   fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
 372   fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
 373   fListHist->Add(fHistDiMuonMass);
 374
 375   fHistDiLeptonMass = new TH1D("fHistDiLeptonMass","Invariant mass of J/#psi candidates",130,2.1,6.0);
 376   fHistDiLeptonMass->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}) (GeV/c)");
 377   fListHist->Add(fHistDiLeptonMass);
 378
 379   TString CutNamePsi2s[14] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Four good tracks",
 380                                 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
 381
 382   fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",14,0.5,14.5);
 383   for (Int_t i = 0; i<14; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
 384   fListHist->Add(fHistNeventsPsi2s);
 385
 386   fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
 387   fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
 388   fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 389   fListHist->Add(fHistPsi2sMassVsPt);
 390
 391   fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",50,2.5,5.5);
 392   fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
 393   fListHist->Add(fHistPsi2sMassCoherent);
 394
 395   fListSystematics = new TList();
 396   fListSystematics->SetOwner();
 397   fListSystematics->SetName("fListSystematics");
 398   fListHist->Add(fListSystematics);
 399   InitSystematics();
 400
 401
 402   PostData(1, fJPsiTree);
 403   PostData(2, fPsi2sTree);
 404   PostData(3, fListTrig);
 405   PostData(4, fListHist);
 406
 407 }//UserCreateOutputObjects
 408
 409 //_____________________________________________________________________________
 410 void AliAnalysisTaskUpcPsi2s::InitSystematics()
 411 {
 412
 413 fListJPsiLoose = new TList();
 414 fListJPsiLoose->SetOwner();
 415 fListJPsiLoose->SetName("JPsiLoose");
 416 fListSystematics->Add(fListJPsiLoose);
 417
 418 TH1D *fHistJPsiNClusLoose = new TH1D("JPsiNClusLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 419 fListJPsiLoose->Add(fHistJPsiNClusLoose);
 420
 421 TH1D *fHistJPsiChi2Loose = new TH1D("JPsiChi2Loose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 422 fListJPsiLoose->Add(fHistJPsiChi2Loose);
 423
 424 TH1D *fHistJPsiDCAzLoose = new TH1D("JPsiDCAzLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 425 fListJPsiLoose->Add(fHistJPsiDCAzLoose);
 426
 427 TH1D *fHistJPsiDCAxyLoose = new TH1D("JPsiDCAxyLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 428 fListJPsiLoose->Add(fHistJPsiDCAxyLoose);
 429
 430 TH1D *fHistJPsiITShitsLoose = new TH1D("JPsiITShitsLoose","Invariant mass of J/#psi candidates",130,2.1,6.0);
 431 fListJPsiLoose->Add(fHistJPsiITShitsLoose);
 432
 433
 434 fListJPsiTight = new TList();
 435 fListJPsiTight->SetOwner();
 436 fListJPsiTight->SetName("JPsiTight");
 437 fListSystematics->Add(fListJPsiTight);
 438
 439 TH1D *fHistJPsiNClusTight = new TH1D("JPsiNClusTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
 440 fListJPsiTight->Add(fHistJPsiNClusTight);
 441
 442 TH1D *fHistJPsiChi2Tight = new TH1D("JPsiChi2Tight","Invariant mass of J/#psi candidates",130,2.1,6.0);
 443 fListJPsiTight->Add(fHistJPsiChi2Tight);
 444
 445 TH1D *fHistJPsiDCAzTight = new TH1D("JPsiDCAzTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
 446 fListJPsiTight->Add(fHistJPsiDCAzTight);
 447
 448 TH1D *fHistJPsiDCAxyTight = new TH1D("JPsiDCAxyTight","Invariant mass of J/#psi candidates",130,2.1,6.0);
 449 fListJPsiTight->Add(fHistJPsiDCAxyTight);
 450
 451
 452 fListPsi2sLoose = new TList();
 453 fListPsi2sLoose->SetOwner();
 454 fListPsi2sLoose->SetName("Psi2sLoose");
 455 fListSystematics->Add(fListPsi2sLoose);
 456
 457 TH1D *fHistPsi2sNClusLoose = new TH1D("Psi2sNClusLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 458 fListPsi2sLoose->Add(fHistPsi2sNClusLoose);
 459
 460 TH1D *fHistPsi2sChi2Loose = new TH1D("Psi2sChi2Loose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 461 fListPsi2sLoose->Add(fHistPsi2sChi2Loose);
 462
 463 TH1D *fHistPsi2sDCAzLoose = new TH1D("Psi2sDCAzLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 464 fListPsi2sLoose->Add(fHistPsi2sDCAzLoose);
 465
 466 TH1D *fHistPsi2sDCAxyLoose = new TH1D("Psi2sDCAxyLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 467 fListPsi2sLoose->Add(fHistPsi2sDCAxyLoose);
 468
 469 TH1D *fHistPsi2sITShitsLoose = new TH1D("Psi2sITShitsLoose","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 470 fListPsi2sLoose->Add(fHistPsi2sITShitsLoose);
 471
 472
 473 fListPsi2sTight = new TList();
 474 fListPsi2sTight->SetOwner();
 475 fListPsi2sTight->SetName("Psi2sTight");
 476 fListSystematics->Add(fListPsi2sTight);
 477
 478 TH1D *fHistPsi2sNClusTight = new TH1D("Psi2sNClusTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 479 fListPsi2sTight->Add(fHistPsi2sNClusTight);
 480
 481 TH1D *fHistPsi2sChi2Tight = new TH1D("Psi2sChi2Tight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 482 fListPsi2sTight->Add(fHistPsi2sChi2Tight);
 483
 484 TH1D *fHistPsi2sDCAzTight = new TH1D("Psi2sDCAzTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 485 fListPsi2sTight->Add(fHistPsi2sDCAzTight);
 486
 487 TH1D *fHistPsi2sDCAxyTight = new TH1D("Psi2sDCAxyTight","Invariant mass of #psi(2S) candidates",50,2.5,5.5);
 488 fListPsi2sTight->Add(fHistPsi2sDCAxyTight);
 489
 490
 491 }
 492
 493 //_____________________________________________________________________________
 494 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
 495 {
 496
 497   //cout<<"#################### Next event ##################"<<endl;
 498
 499   if( fType == 0 ){
 500         RunESDtrig();
 501         if(fRunHist) RunESDhist();
 502         if(fRunTree) RunESDtree();
 503         }
 504
 505   if( fType == 1 ){
 506         RunAODtrig();
 507         if(fRunHist) RunAODhist();
 508         if(fRunTree) RunAODtree();
 509         }
 510
 511 }//UserExec
 512 //_____________________________________________________________________________
 513 void AliAnalysisTaskUpcPsi2s::RunAODtrig()
 514 {
 515
 516   //input event
 517   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 518   if(!aod) return;
 519
 520   fRunNum = aod ->GetRunNumber();
 521   //Trigger
 522   TString trigger = aod->GetFiredTriggerClasses();
 523
 524   if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
 525   if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
 526   if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
 527
 528   if(trigger.Contains("CINT1")) fHistCint1TriggersPerRun->Fill(fRunNum); //CINT1 triggers
 529
 530   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
 531   if(trigger.Contains("CINT1") && (fL0inputs & (1 << 3))) fHistC0tvxAndCint1TriggersPerRun->Fill(fRunNum); //0TVX triggers in CINT1 events
 532
 533   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
 534   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
 535
 536   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 537   if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
 538
 539   //MB, Central and SemiCentral triggers
 540   AliCentrality *centrality = aod->GetCentrality();
 541   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 542
 543   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
 544   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
 545
 546   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
 547
 548   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
 549
 550   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
 551
 552 PostData(3, fListTrig);
 553
 554 }
 555 //_____________________________________________________________________________
 556 void AliAnalysisTaskUpcPsi2s::RunAODhist()
 557 {
 558
 559   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 560
 561   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 562   Double_t muonMass = partMuon->Mass();
 563
 564   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 565   Double_t electronMass = partElectron->Mass();
 566
 567   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 568   Double_t pionMass = partPion->Mass();
 569
 570   //input event
 571   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 572   if(!aod) return;
 573
 574  // cout<<"Event number: "<<((TTree*) GetInputData(0))->GetTree()->GetReadEntry()<<endl;
 575
 576   fHistNeventsJPsi->Fill(1);
 577   fHistNeventsPsi2s->Fill(1);
 578
 579   //Trigger
 580   TString trigger = aod->GetFiredTriggerClasses();
 581
 582   if(!isMC && !trigger.Contains("CCUP") ) return;
 583
 584   fHistNeventsJPsi->Fill(2);
 585   fHistNeventsPsi2s->Fill(2);
 586
 587   //primary vertex
 588   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 589   fVtxContrib = fAODVertex->GetNContributors();
 590   if(fVtxContrib < 2) return;
 591
 592   fHistNeventsJPsi->Fill(3);
 593   fHistNeventsPsi2s->Fill(3);
 594
 595   //VZERO, ZDC
 596   AliAODVZERO *fV0data = aod ->GetVZEROData();
 597   AliAODZDC *fZDCdata = aod->GetZDCData();
 598
 599   fV0Adecision = fV0data->GetV0ADecision();
 600   fV0Cdecision = fV0data->GetV0CDecision();
 601   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
 602
 603   fHistNeventsJPsi->Fill(4);
 604   fHistNeventsPsi2s->Fill(4);
 605
 606   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 607   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 608
 609   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
 610
 611   fHistNeventsJPsi->Fill(5);
 612   fHistNeventsPsi2s->Fill(5);
 613
 614   //Systematics - cut variation
 615   if(fRunSystematics) RunAODsystematics(aod);
 616
 617   //Two tracks loop
 618   Int_t nGoodTracks = 0;
 619   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 620
 621   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
 622   Short_t qLepton[4], qPion[4];
 623   UInt_t nLepton=0, nPion=0, nHighPt=0, nSpdHits=0;
 624   Double_t fRecTPCsignal[5], fRecTPCsignalDist;
 625   Int_t fChannel = 0;
 626   Int_t mass[3]={-1,-1,-1};
 627   Double_t TrackPt[5]={0,0,0,0,0};
 628   Double_t MeanPt = -1;
 629
 630
 631   //Four track loop
 632   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 633     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
 634     if( !trk ) continue;
 635     if(!(trk->TestFilterBit(1<<0))) continue;
 636
 637       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 638       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 639       if(trk->GetTPCNcls() < 50)continue;
 640       if(trk->Chi2perNDF() > 4)continue;
 641       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 642       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 643       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 644       delete trk_clone;
 645       if(TMath::Abs(dca[1]) > 2) continue;
 646       Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 647       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
 648       if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
 649
 650       TrackIndex[nGoodTracks] = itr;
 651       TrackPt[nGoodTracks] = trk->Pt();
 652       nGoodTracks++;
 653
 654       if(nGoodTracks > 4) break;
 655   }//Track loop
 656
 657   nLepton=0; nPion=0; nHighPt=0;
 658   mass[0]= -1; mass[1]= -1, mass[2]= -1;
 659
 660   if(nGoodTracks == 4 && nSpdHits>1){
 661           MeanPt = GetMedian(TrackPt);
 662           fHistNeventsPsi2s->Fill(6);
 663           for(Int_t i=0; i<4; i++){
 664                 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[i]));
 665                 if(!trk) AliFatal("Not a standard AOD");
 666
 667                 if(trk->Pt() > MeanPt){
 668                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 669                         qLepton[nLepton] = trk->Charge();
 670                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 671                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 672                                         mass[nLepton] = 0;
 673                                         }
 674                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 675                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 676                                         mass[nLepton] = 1;
 677                                         }
 678                         nLepton++;
 679                         }
 680                 else{
 681                         qPion[nPion] = trk->Charge();
 682                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
 683                         nPion++;
 684                         }
 685
 686                 if(nLepton > 2 || nPion > 2) break;
 687                 }
 688         if((nLepton == 2) && (nPion == 2)){
 689                 fHistNeventsPsi2s->Fill(7);
 690                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(8);
 691                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(9);
 692                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
 693                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
 694                         fHistNeventsPsi2s->Fill(11);
 695                         if(mass[0] != -1 && mass[1] != -1) {
 696                                 fHistNeventsPsi2s->Fill(12);
 697                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
 698                                 vDilepton = vLepton[0]+vLepton[1];
 699                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
 700                                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
 701                                 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
 702                                 else {
 703                                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
 704                                         if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
 705                                         }
 706
 707                                 if(fChannel == -1) {
 708                                         fHistNeventsPsi2s->Fill(13);
 709                                         if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
 710                                         }
 711                                 if(fChannel == 1){
 712                                         fHistNeventsPsi2s->Fill(14);
 713                                         if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) fHistPsi2sMassCoherent->Fill(vCandidate.M());
 714                                         }
 715                                 }
 716                         }
 717                 }
 718   }
 719
 720   nGoodTracks = 0;
 721   //Two track loop
 722   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 723     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
 724     if( !trk ) continue;
 725     if(!(trk->TestFilterBit(1<<0))) continue;
 726
 727       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 728       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 729       if(trk->GetTPCNcls() < 70)continue;
 730       if(trk->Chi2perNDF() > 4)continue;
 731       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 732       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 733       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 734       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 735       delete trk_clone;
 736       if(TMath::Abs(dca[1]) > 2) continue;
 737       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 738       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
 739
 740       TrackIndex[nGoodTracks] = itr;
 741       nGoodTracks++;
 742
 743       if(nGoodTracks > 2) break;
 744   }//Track loop
 745
 746    nLepton=0; nPion=0; nHighPt=0;
 747   mass[0]= -1; mass[1]= -1, mass[2]= -1;
 748
 749   if(nGoodTracks == 2){
 750           fHistNeventsJPsi->Fill(6);
 751           for(Int_t i=0; i<2; i++){
 752                 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[i]));
 753                 if(!trk) AliFatal("Not a standard AOD");
 754                 if(trk->Pt() > 1) nHighPt++;
 755                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 756                 qLepton[nLepton] = trk->Charge();
 757                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 758                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 759                                 mass[nLepton] = 0;
 760                                 }
 761                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 762                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 763                                 mass[nLepton] = 1;
 764                                 }
 765                 nLepton++;
 766                 }
 767         if(nLepton == 2){
 768                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(7);
 769                 if(qLepton[0]*qLepton[1] < 0){
 770                         fHistNeventsJPsi->Fill(8);
 771                         if(nHighPt > 0){
 772                                 fHistNeventsJPsi->Fill(9);
 773                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
 774                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(10);
 775                                 if(mass[0] != -1 && mass[1] != -1) {
 776                                         fHistNeventsJPsi->Fill(11);
 777                                         vCandidate = vLepton[0]+vLepton[1];
 778                                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
 779                                         if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
 780                                         else {
 781                                                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
 782                                                 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
 783                                                 }
 784                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
 785                                         if(fChannel == -1) {
 786                                                 fHistDiMuonMass->Fill(vCandidate.M());
 787                                                 if(vCandidate.Pt()<0.15)fHistDiLeptonMass->Fill(vCandidate.M());
 788                                                 fHistNeventsJPsi->Fill(12);
 789                                                 }
 790                                         if(fChannel == 1) {
 791                                                 fHistDiElectronMass->Fill(vCandidate.M());
 792                                                 if(vCandidate.Pt()<0.3)fHistDiLeptonMass->Fill(vCandidate.M());
 793                                                 fHistNeventsJPsi->Fill(13);
 794                                                 }
 795                                         }
 796                                 }
 797                         }
 798                 }
 799   }
 800
 801
 802   PostData(4, fListHist);
 803
 804 }
 805
 806 //_____________________________________________________________________________
 807 void AliAnalysisTaskUpcPsi2s::RunAODtree()
 808 {
 809   //input event
 810   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 811   if(!aod) return;
 812
 813   if(isMC) RunAODMC(aod);
 814
 815   //input data
 816   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
 817   fDataFilnam->Clear();
 818   fDataFilnam->SetString(filnam);
 819   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
 820   fRunNum = aod ->GetRunNumber();
 821
 822   //Trigger
 823   TString trigger = aod->GetFiredTriggerClasses();
 824
 825   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
 826   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
 827   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
 828   fTrigger[3]  = trigger.Contains("CINT1");
 829
 830   Bool_t isTriggered = kFALSE;
 831   for(Int_t i=0; i<ntrg; i++) {
 832     if( fTrigger[i] ) isTriggered = kTRUE;
 833   }
 834   if(!isMC && !isTriggered ) return;
 835
 836   //trigger inputs
 837   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
 838   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 839
 840   //Event identification
 841   fPerNum = aod ->GetPeriodNumber();
 842   fOrbNum = aod ->GetOrbitNumber();
 843   fBCrossNum = aod ->GetBunchCrossNumber();
 844
 845   //primary vertex
 846   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 847   fVtxContrib = fAODVertex->GetNContributors();
 848   fVtxPos[0] = fAODVertex->GetX();
 849   fVtxPos[1] = fAODVertex->GetY();
 850   fVtxPos[2] = fAODVertex->GetZ();
 851   Double_t CovMatx[6];
 852   fAODVertex->GetCovarianceMatrix(CovMatx);
 853   fVtxErr[0] = CovMatx[0];
 854   fVtxErr[1] = CovMatx[1];
 855   fVtxErr[2] = CovMatx[2];
 856   fVtxChi2 = fAODVertex->GetChi2();
 857   fVtxNDF = fAODVertex->GetNDF();
 858
 859   //SPD primary vertex
 860   AliAODVertex *fSPDVertex = aod->GetPrimaryVertexSPD();
 861   fSpdVtxContrib = fSPDVertex->GetNContributors();
 862   fSpdVtxPos[0] = fSPDVertex->GetX();
 863   fSpdVtxPos[1] = fSPDVertex->GetY();
 864   fSpdVtxPos[2] = fSPDVertex->GetZ();
 865
 866   //Tracklets
 867   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
 868
 869   //VZERO, ZDC
 870   AliAODVZERO *fV0data = aod ->GetVZEROData();
 871   AliAODZDC *fZDCdata = aod->GetZDCData();
 872
 873   fV0Adecision = fV0data->GetV0ADecision();
 874   fV0Cdecision = fV0data->GetV0CDecision();
 875   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 876   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 877
 878   fZDCAtime = fZDCdata->GetZNATime();
 879   fZDCCtime = fZDCdata->GetZNCTime();
 880
 881   fNLooseTracks = 0;
 882
 883   //Track loop - loose cuts
 884   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 885     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
 886     if( !trk ) continue;
 887     if(!(trk->TestFilterBit(1<<0))) continue;
 888
 889       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 890       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 891       if(trk->GetTPCNcls() < 20)continue;
 892       fNLooseTracks++;
 893   }//Track loop -loose cuts
 894
 895   Int_t nGoodTracks=0;
 896   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 897
 898   //Two track loop
 899   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 900     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
 901     if( !trk ) continue;
 902     if(!(trk->TestFilterBit(1<<0))) continue;
 903
 904       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 905       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 906       if(trk->GetTPCNcls() < 70)continue;
 907       if(trk->Chi2perNDF() > 4)continue;
 908       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 909       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 910       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 911       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 912       delete trk_clone;
 913       if(TMath::Abs(dca[1]) > 2) continue;
 914       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 915       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
 916
 917
 918       TrackIndex[nGoodTracks] = itr;
 919       nGoodTracks++;
 920
 921       if(nGoodTracks > 2) break;
 922   }//Track loop
 923
 924   fJPsiAODTracks->Clear("C");
 925   if(nGoodTracks == 2){
 926
 927           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 928           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 929           Double_t muonMass = partMuon->Mass();
 930           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 931           Double_t electronMass = partElectron->Mass();
 932           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 933           Double_t pionMass = partPion->Mass();
 934
 935           Double_t KFcov[21];
 936           Double_t KFpar[6];
 937           Double_t KFmass = pionMass;
 938           Double_t fRecTPCsignal;
 939           AliKFParticle *KFpart[2];
 940           AliKFVertex *KFvtx = new AliKFVertex();
 941           KFvtx->SetField(aod->GetMagneticField());
 942
 943           for(Int_t i=0; i<2; i++){
 944                 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[i]));
 945                 if(!trk) AliFatal("Not a standard AOD");
 946
 947                 if(fAODVertex->HasDaughter(trk) && trk->GetUsedForVtxFit())fIsVtxContributor[i] = kTRUE;
 948                 else fIsVtxContributor[i] = kFALSE;
 949
 950                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 951                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 952                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 953                 delete trk_clone;
 954
 955                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
 956                 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
 957
 958                 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
 959                 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
 960                 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
 961                 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
 962                 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
 963
 964                 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
 965                 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
 966                 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
 967                 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
 968                 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
 969
 970                 trk->GetPosition(KFpar);
 971                 trk->PxPyPz(KFpar+3);
 972                 trk->GetCovMatrix(KFcov);
 973
 974                 if(trk->Pt() > 1){
 975                         fRecTPCsignal = trk->GetTPCsignal();
 976                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
 977                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
 978                         }
 979                 else KFmass = pionMass;
 980
 981                 KFpart[i] = new AliKFParticle();
 982                 KFpart[i]->SetField(aod->GetMagneticField());
 983                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
 984                 KFvtx->AddDaughter(*KFpart[i]);
 985
 986
 987                 Double_t pos[3]={0,0,0};
 988                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
 989                 parTrk->CopyFromVTrack((AliVTrack*) trk);
 990                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
 991                 else {
 992                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
 993                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
 994                      }
 995                 delete parTrk;
 996                 }
 997   fKfVtxPos[0]= KFvtx->GetX();
 998   fKfVtxPos[1]= KFvtx->GetY();
 999   fKfVtxPos[2]= KFvtx->GetZ();
1000   for(UInt_t i=0; i<2; i++)delete KFpart[i];
1001   delete KFvtx;
1002
1003   if(!isMC) fJPsiTree ->Fill();
1004   }
1005
1006    nGoodTracks = 0;
1007    //Four track loop
1008   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1009     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
1010     if( !trk ) continue;
1011     if(!(trk->TestFilterBit(1<<0))) continue;
1012
1013       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
1014       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
1015       if(trk->GetTPCNcls() < 50)continue;
1016       if(trk->Chi2perNDF() > 4)continue;
1017       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1018       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1019       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1020       delete trk_clone;
1021       if(TMath::Abs(dca[1]) > 2) continue;
1022       Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1023       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
1024
1025       TrackIndex[nGoodTracks] = itr;
1026       nGoodTracks++;
1027
1028       if(nGoodTracks > 4) break;
1029   }//Track loop
1030
1031   fPsi2sAODTracks->Clear("C");
1032   if(nGoodTracks == 4){
1033
1034           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1035           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1036           Double_t muonMass = partMuon->Mass();
1037           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1038           Double_t electronMass = partElectron->Mass();
1039           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1040           Double_t pionMass = partPion->Mass();
1041
1042           Double_t KFcov[21];
1043           Double_t KFpar[6];
1044           Double_t KFmass = pionMass;
1045           Double_t fRecTPCsignal;
1046           AliKFParticle *KFpart[4];
1047           AliKFVertex *KFvtx = new AliKFVertex();
1048           KFvtx->SetField(aod->GetMagneticField());
1049
1050           for(Int_t i=0; i<4; i++){
1051                 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[i]));
1052                 if(!trk) AliFatal("Not a standard AOD");
1053
1054                 if(fAODVertex->HasDaughter(trk) && trk->GetUsedForVtxFit())fIsVtxContributor[i] = kTRUE;
1055                 else fIsVtxContributor[i] = kFALSE;
1056
1057                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1058                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1059                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1060                 delete trk_clone;
1061
1062                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
1063                 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
1064
1065
1066                 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1067                 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1068                 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1069                 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1070                 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
1071
1072                 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1073                 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1074                 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1075                 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1076                 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
1077
1078                 trk->GetPosition(KFpar);
1079                 trk->PxPyPz(KFpar+3);
1080                 trk->GetCovMatrix(KFcov);
1081
1082                 if(trk->Pt() > 1){
1083                         fRecTPCsignal = trk->GetTPCsignal();
1084                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
1085                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
1086                         }
1087                 else KFmass = pionMass;
1088
1089                 KFpart[i] = new AliKFParticle();
1090                 KFpart[i]->SetField(aod->GetMagneticField());
1091                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
1092                 KFvtx->AddDaughter(*KFpart[i]);
1093
1094                 Double_t pos[3]={0,0,0};
1095                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
1096                 parTrk->CopyFromVTrack((AliVTrack*) trk);
1097                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
1098                 else {
1099                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1100                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1101                      }
1102                 delete parTrk;
1103                 }
1104   fKfVtxPos[0]= KFvtx->GetX();
1105   fKfVtxPos[1]= KFvtx->GetY();
1106   fKfVtxPos[2]= KFvtx->GetZ();
1107   for(UInt_t i=0; i<4; i++)delete KFpart[i];
1108   delete KFvtx;
1109   if(!isMC) fPsi2sTree ->Fill();
1110   }
1111
1112   if(isMC){
1113         fJPsiTree ->Fill();
1114         fPsi2sTree ->Fill();
1115   }
1116
1117   PostData(1, fJPsiTree);
1118   PostData(2, fPsi2sTree);
1119
1120 }//RunAOD
1121
1122
1123 //_____________________________________________________________________________
1124 void AliAnalysisTaskUpcPsi2s::RunAODMC(AliAODEvent *aod)
1125 {
1126
1127   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
1128   fTriggerInputsMC[0] = kFALSE;//0SM2
1129   fTriggerInputsMC[1] = fL0inputs & (1 << 0);//0VBA
1130   fTriggerInputsMC[2] = fL0inputs & (1 << 1);//0VBC
1131   fTriggerInputsMC[3] = fL0inputs & (1 << 9);//0OMU
1132
1133   fGenPart->Clear("C");
1134
1135   TClonesArray *arrayMC = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1136   if(!arrayMC) return;
1137
1138   Int_t nmc=0;
1139   //loop over mc particles
1140   for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
1141     AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(imc);
1142     if(!mcPart) continue;
1143
1144     if(mcPart->GetMother() >= 0) continue;
1145
1146     TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1147     part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1148     part->SetPdgCode(mcPart->GetPdgCode());
1149     part->SetUniqueID(imc);
1150   }//loop over mc particles
1151
1152 }//RunAODMC
1153
1154
1155 //_____________________________________________________________________________
1156 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
1157 {
1158
1159   //input event
1160   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1161   if(!esd) return;
1162
1163   fRunNum = esd ->GetRunNumber();
1164   //Trigger
1165   TString trigger = esd->GetFiredTriggerClasses();
1166
1167   if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
1168   if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
1169   if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
1170
1171   if(trigger.Contains("CINT1")) fHistCint1TriggersPerRun->Fill(fRunNum); //CINT1 triggers
1172   if(trigger.Contains("CINT1") && trigger.Contains("C0TVX")) fHistC0tvxAndCint1TriggersPerRun->Fill(fRunNum); //C0TVX triggers in CINT1 events
1173
1174   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
1175   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
1176
1177   if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
1178
1179    //MB, Central and SemiCentral triggers
1180   AliCentrality *centrality = esd->GetCentrality();
1181   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1182
1183   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
1184   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
1185
1186   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
1187
1188   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
1189
1190   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
1191
1192
1193 PostData(3, fListTrig);
1194
1195 }
1196 //_____________________________________________________________________________
1197 void AliAnalysisTaskUpcPsi2s::RunESDhist()
1198 {
1199
1200   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1201
1202   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1203   Double_t muonMass = partMuon->Mass();
1204
1205   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1206   Double_t electronMass = partElectron->Mass();
1207
1208   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1209   Double_t pionMass = partPion->Mass();
1210
1211   //input event
1212   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1213   if(!esd) return;
1214
1215   fHistNeventsJPsi->Fill(1);
1216   fHistNeventsPsi2s->Fill(1);
1217
1218   //Trigger
1219   TString trigger = esd->GetFiredTriggerClasses();
1220
1221   if(!isMC && !trigger.Contains("CCUP") ) return;
1222
1223   fHistNeventsJPsi->Fill(2);
1224   fHistNeventsPsi2s->Fill(2);
1225
1226   //primary vertex
1227   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1228   fVtxContrib = fESDVertex->GetNContributors();
1229   if(fVtxContrib < 2) return;
1230
1231   fHistNeventsJPsi->Fill(3);
1232   fHistNeventsPsi2s->Fill(3);
1233
1234   //VZERO, ZDC
1235   AliESDVZERO *fV0data = esd->GetVZEROData();
1236   AliESDZDC *fZDCdata = esd->GetESDZDC();
1237
1238   fV0Adecision = fV0data->GetV0ADecision();
1239   fV0Cdecision = fV0data->GetV0CDecision();
1240   if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
1241
1242   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1243   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1244   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
1245
1246   fHistNeventsJPsi->Fill(4);
1247   fHistNeventsPsi2s->Fill(4);
1248
1249    Int_t nGoodTracks=0;
1250   //Two tracks loop
1251   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1252
1253   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1254   Short_t qLepton[4], qPion[4];
1255   UInt_t nLepton=0, nPion=0, nHighPt=0;
1256   Double_t fRecTPCsignal[5];
1257   Int_t mass[3]={-1,-1,-1};
1258
1259  //Two Track loop
1260   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1261     AliESDtrack *trk = esd->GetTrack(itr);
1262     if( !trk ) continue;
1263
1264       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1265       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1266       if(trk->GetTPCNcls() < 70)continue;
1267       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1268       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1269       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1270       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1271       trk->GetImpactParameters(dca[0],dca[1]);
1272       if(TMath::Abs(dca[1]) > 2) continue;
1273       if(TMath::Abs(dca[1]) > 0.2) continue;
1274
1275       TrackIndex[nGoodTracks] = itr;
1276       nGoodTracks++;
1277       if(nGoodTracks > 2) break;
1278   }//Track loop
1279
1280
1281   if(nGoodTracks == 2){
1282           fHistNeventsJPsi->Fill(5);
1283           for(Int_t i=0; i<2; i++){
1284                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1285                 if(trk->Pt() > 1) nHighPt++;
1286                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1287                 qLepton[nLepton] = trk->Charge();
1288                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1289                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1290                                 mass[nLepton] = 0;
1291                                 }
1292                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1293                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1294                                 mass[nLepton] = 1;
1295                                 }
1296                 nLepton++;
1297                 }
1298         if(nLepton == 2){
1299                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
1300                 if(qLepton[0]*qLepton[1] < 0){
1301                         fHistNeventsJPsi->Fill(7);
1302                         if(nHighPt > 0){
1303                                 fHistNeventsJPsi->Fill(8);
1304                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
1305                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
1306                                 if(mass[0] == mass[1] && mass[0] != -1) {
1307                                         fHistNeventsJPsi->Fill(10);
1308                                         vCandidate = vLepton[0]+vLepton[1];
1309                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
1310                                         if(mass[0] == 0) {
1311                                                 fHistDiMuonMass->Fill(vCandidate.M());
1312                                                 fHistNeventsJPsi->Fill(11);
1313                                                 }
1314                                         if(mass[0] == 1) {
1315                                                 fHistDiElectronMass->Fill(vCandidate.M());
1316                                                 fHistNeventsJPsi->Fill(12);
1317                                                 }
1318                                         }
1319                                 }
1320                         }
1321                 }
1322   }
1323   nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1324   mass[0]= -1; mass[1]= -1, mass[2]= -1;
1325
1326     //Four Track loop
1327   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1328     AliESDtrack *trk = esd->GetTrack(itr);
1329     if( !trk ) continue;
1330
1331       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1332       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1333       if(trk->GetTPCNcls() < 50)continue;
1334       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1335       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1336       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1337       trk->GetImpactParameters(dca[0],dca[1]);
1338       if(TMath::Abs(dca[1]) > 2) continue;
1339
1340       TrackIndex[nGoodTracks] = itr;
1341       nGoodTracks++;
1342       if(nGoodTracks > 4) break;
1343   }//Track loop
1344
1345   if(nGoodTracks == 4){
1346           fHistNeventsPsi2s->Fill(5);
1347           for(Int_t i=0; i<4; i++){
1348                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1349
1350                 if(trk->Pt() > 1){
1351                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1352                         qLepton[nLepton] = trk->Charge();
1353                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1354                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1355                                         mass[nLepton] = 0;
1356                                         }
1357                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1358                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1359                                         mass[nLepton] = 1;
1360                                         }
1361                         nLepton++;
1362                         }
1363                 else{
1364                         qPion[nPion] = trk->Charge();
1365                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1366                         nPion++;
1367                         }
1368
1369                 if(nLepton > 2 || nPion > 2) break;
1370                 }
1371         if((nLepton == 2) && (nPion == 2)){
1372                 fHistNeventsPsi2s->Fill(6);
1373                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1374                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1375                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1376                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1377                         fHistNeventsPsi2s->Fill(10);
1378                         if(mass[0] == mass[1]) {
1379                                 fHistNeventsPsi2s->Fill(11);
1380                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1381                                 vDilepton = vLepton[0]+vLepton[1];
1382                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1383                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1384                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1385                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1386                                 }
1387                         }
1388                 }
1389   }
1390
1391   PostData(4, fListHist);
1392
1393 }
1394
1395 //_____________________________________________________________________________
1396 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1397 {
1398
1399   //input event
1400   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1401   if(!esd) return;
1402
1403   if(isMC) RunESDMC(esd);
1404
1405   //input data
1406   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1407   fDataFilnam->Clear();
1408   fDataFilnam->SetString(filnam);
1409   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1410   fRunNum = esd->GetRunNumber();
1411
1412    //Trigger
1413   TString trigger = esd->GetFiredTriggerClasses();
1414
1415   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1416   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
1417   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1418
1419   Bool_t isTriggered = kFALSE;
1420   for(Int_t i=0; i<ntrg; i++) {
1421     if( fTrigger[i] ) isTriggered = kTRUE;
1422   }
1423   if(!isMC && !isTriggered ) return;
1424
1425   //trigger inputs
1426   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1427   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1428
1429   //TOF trigger info (0OMU)
1430   fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1431   fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1432
1433   //Event identification
1434   fPerNum = esd->GetPeriodNumber();
1435   fOrbNum = esd->GetOrbitNumber();
1436   fBCrossNum = esd->GetBunchCrossNumber();
1437
1438   //primary vertex
1439   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1440   fVtxContrib = fESDVertex->GetNContributors();
1441   fVtxPos[0] = fESDVertex->GetX();
1442   fVtxPos[1] = fESDVertex->GetY();
1443   fVtxPos[2] = fESDVertex->GetZ();
1444   Double_t CovMatx[6];
1445   fESDVertex->GetCovarianceMatrix(CovMatx);
1446   fVtxErr[0] = CovMatx[0];
1447   fVtxErr[1] = CovMatx[1];
1448   fVtxErr[2] = CovMatx[2];
1449   fVtxChi2 = fESDVertex->GetChi2();
1450   fVtxNDF = fESDVertex->GetNDF();
1451
1452   //SPD primary vertex
1453   AliESDVertex *fSPDVertex = (AliESDVertex*) esd->GetPrimaryVertexSPD();
1454   if(fSPDVertex){
1455         fSpdVtxPos[0] = fSPDVertex->GetX();
1456         fSpdVtxPos[1] = fSPDVertex->GetY();
1457         fSpdVtxPos[2] = fSPDVertex->GetZ();
1458         }
1459
1460   //Tracklets
1461   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1462
1463   //VZERO, ZDC
1464   AliESDVZERO *fV0data = esd->GetVZEROData();
1465   AliESDZDC *fZDCdata = esd->GetESDZDC();
1466
1467   fV0Adecision = fV0data->GetV0ADecision();
1468   fV0Cdecision = fV0data->GetV0CDecision();
1469   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1470   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1471
1472   fNLooseTracks = 0;
1473
1474   //Track loop - loose cuts
1475   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1476     AliESDtrack *trk = esd->GetTrack(itr);
1477     if( !trk ) continue;
1478
1479       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1480       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1481       if(trk->GetTPCNcls() < 20)continue;
1482       fNLooseTracks++;
1483   }//Track loop -loose cuts
1484
1485   Int_t nGoodTracks=0;
1486   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1487
1488   //Two Track loop
1489   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1490     AliESDtrack *trk = esd->GetTrack(itr);
1491     if( !trk ) continue;
1492
1493       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1494       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1495       if(trk->GetTPCNcls() < 70)continue;
1496       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1497       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1498       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1499       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1500       trk->GetImpactParameters(dca[0],dca[1]);
1501       if(TMath::Abs(dca[1]) > 2) continue;
1502       if(TMath::Abs(dca[0]) > 0.2) continue;
1503
1504       TrackIndex[nGoodTracks] = itr;
1505       nGoodTracks++;
1506       if(nGoodTracks > 2) break;
1507   }//Track loop
1508
1509   fJPsiESDTracks->Clear("C");
1510   if(nGoodTracks == 2){
1511           for(Int_t i=0; i<2; i++){
1512                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1513
1514                 AliExternalTrackParam cParam;
1515                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1516
1517                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
1518
1519                 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1520                 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1521                 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1522                 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1523                 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
1524
1525                 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1526                 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1527                 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1528                 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1529                 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
1530
1531                 Double_t pos[3]={0,0,0};
1532                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1533                 else {
1534                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1535                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1536                      }
1537                 }
1538   if(!isMC) fJPsiTree ->Fill();
1539   }
1540
1541   nGoodTracks = 0;
1542   //Four track loop
1543   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1544     AliESDtrack *trk = esd->GetTrack(itr);
1545     if( !trk ) continue;
1546
1547       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1548       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1549       if(trk->GetTPCNcls() < 50)continue;
1550       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1551       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1552       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1553       trk->GetImpactParameters(dca[0],dca[1]);
1554       if(TMath::Abs(dca[1]) > 2) continue;
1555
1556       TrackIndex[nGoodTracks] = itr;
1557       nGoodTracks++;
1558       if(nGoodTracks > 4) break;
1559   }//Track loop
1560
1561   fPsi2sESDTracks->Clear("C");
1562   if(nGoodTracks == 4){
1563           for(Int_t i=0; i<4; i++){
1564                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1565
1566                 AliExternalTrackParam cParam;
1567                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1568
1569                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1570
1571                 fPIDTPCMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1572                 fPIDTPCElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1573                 fPIDTPCPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1574                 fPIDTPCKaon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kKaon);
1575                 fPIDTPCProton[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kProton);
1576
1577                 fPIDTOFMuon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kMuon);
1578                 fPIDTOFElectron[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kElectron);
1579                 fPIDTOFPion[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kPion);
1580                 fPIDTOFKaon[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kKaon);
1581                 fPIDTOFProton[i] = fPIDResponse->NumberOfSigmasTOF(trk,AliPID::kProton);
1582
1583                 Double_t pos[3]={0,0,0};
1584                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1585                 else {
1586                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1587                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1588                      }
1589                 }
1590   if(!isMC) fPsi2sTree ->Fill();
1591   }
1592
1593   if(isMC){
1594         fJPsiTree ->Fill();
1595         fPsi2sTree ->Fill();
1596   }
1597
1598   PostData(1, fJPsiTree);
1599   PostData(2, fPsi2sTree);
1600
1601 }//RunESD
1602
1603
1604 //_____________________________________________________________________________
1605 void AliAnalysisTaskUpcPsi2s::RunESDMC(AliESDEvent* esd)
1606 {
1607
1608   AliTriggerAnalysis *fTrigAna = new AliTriggerAnalysis();
1609   fTrigAna->SetAnalyzeMC(isMC);
1610
1611   if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) > 1) fTriggerInputsMC[0] = kTRUE;
1612   if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) < 2) fTriggerInputsMC[0] = kFALSE;
1613
1614   fTriggerInputsMC[1] = esd->GetHeader()->IsTriggerInputFired("0OMU");
1615   fTriggerInputsMC[2] = esd->GetHeader()->IsTriggerInputFired("0VBA");
1616   fTriggerInputsMC[3] = esd->GetHeader()->IsTriggerInputFired("0VBC");
1617
1618   fGenPart->Clear("C");
1619
1620   AliMCEvent *mc = MCEvent();
1621   if(!mc) return;
1622
1623   Int_t nmc = 0;
1624   //loop over mc particles
1625   for(Int_t imc=0; imc<mc->GetNumberOfTracks(); imc++) {
1626     AliMCParticle *mcPart = (AliMCParticle*) mc->GetTrack(imc);
1627     if(!mcPart) continue;
1628
1629     if(mcPart->GetMother() >= 0) continue;
1630
1631     TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1632     part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1633     part->SetPdgCode(mcPart->PdgCode());
1634     part->SetUniqueID(imc);
1635   }//loop over mc particles
1636
1637 }//RunESDMC
1638
1639
1640
1641 //_____________________________________________________________________________
1642 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
1643 {
1644
1645   cout<<"Analysis complete."<<endl;
1646 }//Terminate
1647
1648 //_____________________________________________________________________________
1649 Double_t AliAnalysisTaskUpcPsi2s::GetMedian(Double_t *daArray) {
1650     // Allocate an array of the same size and sort it.
1651     Double_t dpSorted[4];
1652     for (Int_t i = 0; i < 4; ++i) {
1653         dpSorted[i] = daArray[i];
1654     }
1655     for (Int_t i = 3; i > 0; --i) {
1656         for (Int_t j = 0; j < i; ++j) {
1657             if (dpSorted[j] > dpSorted[j+1]) {
1658                 Double_t dTemp = dpSorted[j];
1659                 dpSorted[j] = dpSorted[j+1];
1660                 dpSorted[j+1] = dTemp;
1661             }
1662         }
1663     }
1664
1665     // Middle or average of middle values in the sorted array.
1666     Double_t dMedian = 0.0;
1667     dMedian = (dpSorted[2] + dpSorted[1])/2.0;
1668
1669     return dMedian;
1670 }
1671
1672 //_____________________________________________________________________________
1673 void AliAnalysisTaskUpcPsi2s::RunAODsystematics(AliAODEvent* aod)
1674 {
1675
1676   Double_t fJPsiSels[4];
1677
1678   fJPsiSels[0] =   70; //min number of TPC clusters
1679   fJPsiSels[1] =   4; //chi2
1680   fJPsiSels[2] =   2; //DCAz
1681   fJPsiSels[3] =   1; // DCAxy 1x
1682
1683   Double_t fJPsiSelsMid[4];
1684
1685   fJPsiSelsMid[0] =   70; //min number of TPC clusters
1686   fJPsiSelsMid[1] =   4; //chi2
1687   fJPsiSelsMid[2] =   2; //DCAz
1688   fJPsiSelsMid[3] =   1; // DCAxy 1x
1689
1690   Double_t fJPsiSelsLoose[4];
1691
1692   fJPsiSelsLoose[0] =   60; //min number of TPC clusters
1693   fJPsiSelsLoose[1] =   5; //chi2
1694   fJPsiSelsLoose[2] =   3; //DCAz
1695   fJPsiSelsLoose[3] =   2; // DCAxy 2x
1696
1697   Double_t fJPsiSelsTight[4];
1698
1699   fJPsiSelsTight[0] =   80; //min number of TPC clusters
1700   fJPsiSelsTight[1] =   3.5; //chi2
1701   fJPsiSelsTight[2] =   1; //DCAz
1702   fJPsiSelsTight[3] =   0.5; // DCAxy 0.5x
1703
1704   Int_t nGoodTracks = 0;
1705   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1706
1707   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1708   Short_t qLepton[4],qPion[4];
1709   UInt_t nLepton=0, nPion=0, nHighPt=0;
1710   Double_t fRecTPCsignal[5], fRecTPCsignalDist;
1711   Int_t fChannel = 0;
1712
1713   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
1714
1715   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
1716
1717   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
1718   Double_t muonMass = partMuon->Mass();
1719
1720   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
1721   Double_t electronMass = partElectron->Mass();
1722
1723   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
1724   Double_t pionMass = partPion->Mass();
1725
1726
1727 for(Int_t i=0; i<5; i++){
1728           //cout<<"Loose sytematics, cut"<<i<<endl;
1729           for(Int_t j=0; j<4; j++){
1730                   if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
1731                   else fJPsiSels[j] = fJPsiSelsMid[j];
1732           }
1733   //Two track loop
1734   nGoodTracks = 0;
1735   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1736     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
1737     if( !trk ) continue;
1738     if(!(trk->TestFilterBit(1<<0))) continue;
1739
1740       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1741       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1742       if(i!=4){ if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;}
1743       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1744       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1745       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1746       delete trk_clone;
1747       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1748
1749       if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1750       if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1751       if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1752       if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1753
1754       TrackIndex[nGoodTracks] = itr;
1755       nGoodTracks++;
1756
1757       if(nGoodTracks > 2) break;
1758   }//Track loop
1759
1760   Int_t mass[3]={-1,-1,-1};
1761   fChannel = 0;
1762   nLepton=0; nHighPt=0;
1763
1764   if(nGoodTracks == 2){
1765           for(Int_t k=0; k<2; k++){
1766                 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[k]));
1767                 if(!trk) AliFatal("Not a standard AOD");
1768
1769                 if(trk->Pt() > 1) nHighPt++;
1770                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1771                 qLepton[nLepton] = trk->Charge();
1772                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1773                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1774                                 mass[nLepton] = 0;
1775                                 }
1776                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1777                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1778                                 mass[nLepton] = 1;
1779                                 }
1780                 nLepton++;
1781                 }
1782         if(nLepton == 2){
1783                 if(qLepton[0]*qLepton[1] < 0 && nHighPt > 0 && (mass[0]!=-1 || mass[1]!=-1)){
1784                         vCandidate = vLepton[0]+vLepton[1];
1785                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1786                         if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1787                         else {
1788                                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1789                                 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
1790                                 }
1791                         if(fChannel == -1 && vCandidate.Pt()<0.15) ((TH1D*)(fListJPsiLoose->At(i)))->Fill(vCandidate.M());
1792                         if(fChannel == 1 && vCandidate.Pt()<0.3) ((TH1D*)(fListJPsiLoose->At(i)))->Fill(vCandidate.M());
1793                         }
1794                 }
1795   }
1796 }//loose cuts
1797
1798 for(Int_t i=0; i<4; i++){
1799           //cout<<"Tight sytematics, cut"<<i<<endl;
1800           for(Int_t j=0; j<4; j++){
1801                   if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
1802                   else fJPsiSels[j] = fJPsiSelsMid[j];
1803           }
1804   //Two track loop
1805   nGoodTracks = 0;
1806   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1807     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
1808     if( !trk ) continue;
1809     if(!(trk->TestFilterBit(1<<0))) continue;
1810
1811       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1812       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1813       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1814       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1815       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1816       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1817       delete trk_clone;
1818       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1819
1820       if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1821       if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1822       if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1823       if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1824
1825       TrackIndex[nGoodTracks] = itr;
1826       nGoodTracks++;
1827
1828       if(nGoodTracks > 2) break;
1829   }//Track loop
1830
1831   Int_t mass[3]={-1,-1,-1};
1832   fChannel = 0;
1833   nLepton=0; nHighPt=0;
1834
1835   if(nGoodTracks == 2){
1836           for(Int_t k=0; k<2; k++){
1837                 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[k]));
1838                 if(!trk) AliFatal("Not a standard AOD");
1839
1840                 if(trk->Pt() > 1) nHighPt++;
1841                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1842                 qLepton[nLepton] = trk->Charge();
1843                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1844                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1845                                 mass[nLepton] = 0;
1846                                 }
1847                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1848                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1849                                 mass[nLepton] = 1;
1850                                 }
1851                 nLepton++;
1852                 }
1853         if(nLepton == 2){
1854                 if(qLepton[0]*qLepton[1] < 0 && nHighPt > 0 && (mass[0]!=-1 || mass[1]!=-1)){
1855                         vCandidate = vLepton[0]+vLepton[1];
1856                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1857                         if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1858                         else {
1859                                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1860                                 if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
1861                                 }
1862                         if(fChannel == -1 && vCandidate.Pt()<0.15) ((TH1D*)(fListJPsiTight->At(i)))->Fill(vCandidate.M());
1863                         if(fChannel == 1 && vCandidate.Pt()<0.3) ((TH1D*)(fListJPsiTight->At(i)))->Fill(vCandidate.M());
1864                         }
1865                 }
1866   }
1867 }//tight cuts
1868
1869 //---------------------------------------------Psi2s------------------------------------------------------------------------
1870
1871   Double_t fPsi2sSels[4];
1872
1873   fPsi2sSels[0] =   50; //min number of TPC clusters
1874   fPsi2sSels[1] =   4; //chi2
1875   fPsi2sSels[2] =   2; //DCAz
1876   fPsi2sSels[3] =   4; // DCAxy 1x
1877
1878   Double_t fPsi2sSelsMid[4];
1879
1880   fPsi2sSelsMid[0] =   50; //min number of TPC clusters
1881   fPsi2sSelsMid[1] =   4; //chi2
1882   fPsi2sSelsMid[2] =   2; //DCAz
1883   fPsi2sSelsMid[3] =   4; // DCAxy 1x
1884
1885   Double_t fPsi2sSelsLoose[4];
1886
1887   fPsi2sSelsLoose[0] =   60; //min number of TPC clusters
1888   fPsi2sSelsLoose[1] =   5; //chi2
1889   fPsi2sSelsLoose[2] =   3; //DCAz
1890   fPsi2sSelsLoose[3] =   6; // DCAxy 2x
1891
1892   Double_t fPsi2sSelsTight[4];
1893
1894   fPsi2sSelsTight[0] =   70; //min number of TPC clusters
1895   fPsi2sSelsTight[1] =   3.5; //chi2
1896   fPsi2sSelsTight[2] =   1; //DCAz
1897   fPsi2sSelsTight[3] =   2; // DCAxy 0.5x
1898
1899   nGoodTracks = 0; nLepton=0; nHighPt=0; fChannel = 0;
1900   Int_t nSpdHits = 0;
1901   Double_t TrackPt[5]={0,0,0,0,0};
1902   Double_t MeanPt = -1;
1903
1904 for(Int_t i=0; i<5; i++){
1905           //cout<<"Loose systematics psi2s, cut"<<i<<endl;
1906           for(Int_t j=0; j<4; j++){
1907                   if(i==j) fJPsiSels[j] = fJPsiSelsLoose[i];
1908                   else fJPsiSels[j] = fJPsiSelsMid[j];
1909           }
1910
1911   //Four track loop
1912   nGoodTracks = 0; nSpdHits = 0;
1913   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1914     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
1915     if( !trk ) continue;
1916     if(!(trk->TestFilterBit(1<<0))) continue;
1917
1918       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1919       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1920       if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1921       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
1922       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
1923       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1924       delete trk_clone;
1925       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
1926
1927       if(trk->GetTPCNcls() < fJPsiSels[0])continue;
1928       if(trk->Chi2perNDF() > fJPsiSels[1])continue;
1929       if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
1930       if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
1931       if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
1932
1933       TrackIndex[nGoodTracks] = itr;
1934       TrackPt[nGoodTracks] = trk->Pt();
1935       nGoodTracks++;
1936
1937       if(nGoodTracks > 4) break;
1938   }//Track loop
1939
1940   Int_t mass[3]={-1,-1,-1};
1941   fChannel = 0;
1942   nLepton=0; nPion=0; nHighPt=0;
1943
1944   if(nGoodTracks == 4){
1945           if(i!=4){ if(nSpdHits<2) continue;}
1946           MeanPt = GetMedian(TrackPt);
1947           for(Int_t k=0; k<4; k++){
1948                 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[k]));
1949                 if(!trk) AliFatal("Not a standard AOD");
1950
1951                 if(trk->Pt() > MeanPt){
1952                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1953                         qLepton[nLepton] = trk->Charge();
1954                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1955                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1956                                         mass[nLepton] = 0;
1957                                         }
1958                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1959                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1960                                         mass[nLepton] = 1;
1961                                         }
1962                         nLepton++;
1963                         }
1964                 else{
1965                         qPion[nPion] = trk->Charge();
1966                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1967                         nPion++;
1968                         }
1969                 }
1970         if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0) && mass[0] != -1 && mass[1] != -1){
1971                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1972                 vDilepton = vLepton[0]+vLepton[1];
1973                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
1974                 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
1975                 else {
1976                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
1977                         if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
1978                         }
1979                 if(fChannel == -1) if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) ((TH1D*)(fListPsi2sLoose->At(i)))->Fill(vCandidate.M());
1980                 if(fChannel == 1) if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) ((TH1D*)(fListPsi2sLoose->At(i)))->Fill(vCandidate.M());
1981         }
1982   }
1983 }//loose cuts
1984
1985 for(Int_t i=0; i<4; i++){
1986           //cout<<"Tight systematics psi2s, cut"<<i<<endl;
1987           for(Int_t j=0; j<4; j++){
1988                   if(i==j) fJPsiSels[j] = fJPsiSelsTight[i];
1989                   else fJPsiSels[j] = fJPsiSelsMid[j];
1990           }
1991
1992   //Four track loop
1993   nGoodTracks = 0; nSpdHits = 0;
1994   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
1995     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
1996     if( !trk ) continue;
1997     if(!(trk->TestFilterBit(1<<0))) continue;
1998
1999       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
2000       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
2001       if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
2002       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
2003       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
2004       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
2005       delete trk_clone;
2006       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
2007
2008       if(trk->GetTPCNcls() < fJPsiSels[0])continue;
2009       if(trk->Chi2perNDF() > fJPsiSels[1])continue;
2010       if(TMath::Abs(dca[1]) > fJPsiSels[2]) continue;
2011       if(TMath::Abs(dca[0]) > fJPsiSels[3]*cut_DCAxy) continue;
2012       if((trk->HasPointOnITSLayer(0))||(trk->HasPointOnITSLayer(1))) nSpdHits++;
2013
2014       TrackIndex[nGoodTracks] = itr;
2015       TrackPt[nGoodTracks] = trk->Pt();
2016       nGoodTracks++;
2017
2018       if(nGoodTracks > 4) break;
2019   }//Track loop
2020
2021   Int_t mass[3]={-1,-1,-1};
2022   fChannel = 0;
2023     nLepton=0; nPion=0; nHighPt=0;
2024
2025   if(nGoodTracks == 4){
2026           if(nSpdHits<2) continue;
2027           MeanPt = GetMedian(TrackPt);
2028           for(Int_t k=0; k<4; k++){
2029                 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(TrackIndex[k]));
2030                 if(!trk) AliFatal("Not a standard AOD");
2031
2032                 if(trk->Pt() > MeanPt){
2033                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
2034                         qLepton[nLepton] = trk->Charge();
2035                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
2036                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
2037                                         mass[nLepton] = 0;
2038                                         }
2039                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
2040                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
2041                                         mass[nLepton] = 1;
2042                                         }
2043                         nLepton++;
2044                         }
2045                 else{
2046                         qPion[nPion] = trk->Charge();
2047                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
2048                         nPion++;
2049                         }
2050                 }
2051         if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0) && mass[0] != -1 && mass[1] != -1){
2052                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
2053                 vDilepton = vLepton[0]+vLepton[1];
2054                 fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-56,2)+TMath::Power(fRecTPCsignal[1]-56,2));
2055                 if (fRecTPCsignalDist < 3.6*4.0) fChannel = -1;
2056                 else {
2057                         fRecTPCsignalDist = TMath::Sqrt(TMath::Power(fRecTPCsignal[0]-78,2)+TMath::Power(fRecTPCsignal[1]-78,2));
2058                         if (fRecTPCsignalDist < 4.1*4.0) fChannel = 1;
2059                         }
2060                 if(fChannel == -1) if(vDilepton.M() > 3.0 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.15) ((TH1D*)(fListPsi2sTight->At(i)))->Fill(vCandidate.M());
2061                 if(fChannel == 1) if(vDilepton.M() > 2.6 && vDilepton.M() < 3.2 && vCandidate.Pt()<0.3) ((TH1D*)(fListPsi2sTight->At(i)))->Fill(vCandidate.M());
2062         }
2063   }
2064 }//Tight cuts
2065
2066 }