PWGUD/UPC/AliAnalysisTaskUpcPsi2s.cxx

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