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