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