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),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
  69     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
  70     fTOFtrig1(0), fTOFtrig2(0),
  71     fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
  72     fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
  73     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
  74     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
  75     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
  76     fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),
  77     fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
  78     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
  79     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
  80
  81 {
  82
  83 //Dummy constructor
  84
  85 }//AliAnalysisTaskUpcPsi2s
  86
  87
  88 //_____________________________________________________________________________
  89 AliAnalysisTaskUpcPsi2s::AliAnalysisTaskUpcPsi2s(const char *name)
  90   : AliAnalysisTaskSE(name),fType(0),isMC(kFALSE),fRunTree(kTRUE),fRunHist(kTRUE),fPIDResponse(0),fJPsiTree(0),fPsi2sTree(0),
  91     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),
  92     fTOFtrig1(0), fTOFtrig2(0),
  93     fVtxContrib(0),fVtxChi2(0),fVtxNDF(0),
  94     fBCrossNum(0),fNtracklets(0),fNLooseTracks(0),
  95     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
  96     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
  97     fJPsiAODTracks(0),fJPsiESDTracks(0),fPsi2sAODTracks(0),fPsi2sESDTracks(0),fGenPart(0),
  98     fListTrig(0),fHistCcup4TriggersPerRun(0), fHistCcup7TriggersPerRun(0), fHistCcup2TriggersPerRun(0),
  99     fHistZedTriggersPerRun(0),fHistCvlnTriggersPerRun(0), fHistMBTriggersPerRun(0),fHistCentralTriggersPerRun(0),fHistSemiCentralTriggersPerRun(0),
 100     fListHist(0),fHistNeventsJPsi(0),fHistTPCsignalJPsi(0),fHistDiLeptonPtJPsi(0),fHistDiElectronMass(0),fHistDiMuonMass(0),
 101     fHistNeventsPsi2s(0),fHistPsi2sMassVsPt(0),fHistPsi2sMassCoherent(0)
 102
 103 {
 104
 105   // Constructor
 106   if( strstr(name,"ESD") ) fType = 0;
 107   if( strstr(name,"AOD") ) fType = 1;
 108
 109   Init();
 110
 111   DefineOutput(1, TTree::Class());
 112   DefineOutput(2, TTree::Class());
 113   DefineOutput(3, TList::Class());
 114   DefineOutput(4, TList::Class());
 115
 116 }//AliAnalysisTaskUpcPsi2s
 117
 118 //_____________________________________________________________________________
 119 void AliAnalysisTaskUpcPsi2s::Init()
 120 {
 121
 122   for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
 123   for(Int_t i=0; i<4; i++) {
 124         fTOFphi[i] = -666;
 125         fPIDMuon[i] = -666;
 126         fPIDElectron[i] = -666;
 127         fPIDPion[i] = -666;
 128         fTriggerInputsMC[i] = kFALSE;
 129         }
 130   for(Int_t i=0; i<3; i++){
 131         fVtxPos[i] = -666;
 132         fVtxErr[i] = -666;
 133         fKfVtxPos[i] = -666;
 134         }
 135
 136 }//Init
 137
 138 //_____________________________________________________________________________
 139 AliAnalysisTaskUpcPsi2s::~AliAnalysisTaskUpcPsi2s()
 140 {
 141   // Destructor
 142   if(fJPsiTree){
 143      delete fJPsiTree;
 144      fJPsiTree = 0x0;
 145   }
 146   if(fPsi2sTree){
 147      delete fPsi2sTree;
 148      fPsi2sTree = 0x0;
 149   }
 150   if(fListTrig){
 151      delete fListTrig;
 152      fListTrig = 0x0;
 153   }
 154   if(fListHist){
 155      delete fListHist;
 156      fListHist = 0x0;
 157   }
 158
 159 }//~AliAnalysisTaskUpcPsi2s
 160
 161
 162 //_____________________________________________________________________________
 163 void AliAnalysisTaskUpcPsi2s::UserCreateOutputObjects()
 164 {
 165   //PID response
 166   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 167   AliInputEventHandler *inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 168   fPIDResponse = inputHandler->GetPIDResponse();
 169
 170   //input file
 171   fDataFilnam = new TObjString();
 172   fDataFilnam->SetString("");
 173
 174     //tracks
 175   fJPsiAODTracks = new TClonesArray("AliAODTrack", 1000);
 176   fJPsiESDTracks = new TClonesArray("AliESDtrack", 1000);
 177   fPsi2sAODTracks = new TClonesArray("AliAODTrack", 1000);
 178   fPsi2sESDTracks = new TClonesArray("AliESDtrack", 1000);
 179   fGenPart = new TClonesArray("TParticle", 1000);
 180
 181   //output tree with JPsi candidate events
 182   fJPsiTree = new TTree("fJPsiTree", "fJPsiTree");
 183   fJPsiTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 184   fJPsiTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 185   fJPsiTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 186
 187   fJPsiTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 188   fJPsiTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 189   fJPsiTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 190   fJPsiTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 191   fJPsiTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 192   fJPsiTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
 193   fJPsiTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 194
 195   fJPsiTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
 196   fJPsiTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
 197   fJPsiTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
 198
 199   fJPsiTree ->Branch("fPIDMuon", &fPIDMuon[0], "fPIDMuon[4]/D");
 200   fJPsiTree ->Branch("fPIDElectron", &fPIDElectron[0], "fPIDElectron[4]/D");
 201   fJPsiTree ->Branch("fPIDPion", &fPIDPion[0], "fPIDPion[4]/D");
 202
 203   fJPsiTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
 204   fJPsiTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
 205   fJPsiTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
 206   fJPsiTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
 207
 208   fJPsiTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
 209
 210   fJPsiTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 211   fJPsiTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 212   fJPsiTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 213   fJPsiTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 214   fJPsiTree ->Branch("fDataFilnam", &fDataFilnam);
 215   fJPsiTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 216   fJPsiTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 217   if( fType == 0 ) {
 218     fJPsiTree ->Branch("fJPsiESDTracks", &fJPsiESDTracks);
 219   }
 220   if( fType == 1 ) {
 221     fJPsiTree ->Branch("fJPsiAODTracks", &fJPsiAODTracks);
 222   }
 223   if(isMC) {
 224     fJPsiTree ->Branch("fGenPart", &fGenPart);
 225     fJPsiTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
 226   }
 227
 228
 229  //output tree with Psi2s candidate events
 230   fPsi2sTree = new TTree("fPsi2sTree", "fPsi2sTree");
 231   fPsi2sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 232   fPsi2sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 233   fPsi2sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 234
 235   fPsi2sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 236   fPsi2sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 237   fPsi2sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 238   fPsi2sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 239   fPsi2sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 240   fPsi2sTree ->Branch("fNLooseTracks", &fNLooseTracks, "fNLooseTracks/s");
 241   fPsi2sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 242
 243   fPsi2sTree ->Branch("fTOFtrig1", &fTOFtrig1, "fTOFtrig1/O");
 244   fPsi2sTree ->Branch("fTOFtrig2", &fTOFtrig2, "fTOFtrig2/O");
 245   fPsi2sTree ->Branch("fTOFphi", &fTOFphi[0], "fTOFphi[4]/D");
 246
 247   fPsi2sTree ->Branch("fPIDMuon", &fPIDMuon[0], "fPIDMuon[4]/D");
 248   fPsi2sTree ->Branch("fPIDElectron", &fPIDElectron[0], "fPIDElectron[4]/D");
 249   fPsi2sTree ->Branch("fPIDPion", &fPIDPion[0], "fPIDPion[4]/D");
 250
 251   fPsi2sTree ->Branch("fVtxPos", &fVtxPos[0], "fVtxPos[3]/D");
 252   fPsi2sTree ->Branch("fVtxErr", &fVtxErr[0], "fVtxErr[3]/D");
 253   fPsi2sTree ->Branch("fVtxChi2", &fVtxChi2, "fVtxChi2/D");
 254   fPsi2sTree ->Branch("fVtxNDF", &fVtxNDF, "fVtxNDF/D");
 255
 256   fPsi2sTree ->Branch("fKfVtxPos", &fKfVtxPos[0], "fKfVtxPos[3]/D");
 257
 258   fPsi2sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 259   fPsi2sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 260   fPsi2sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 261   fPsi2sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 262   fPsi2sTree ->Branch("fDataFilnam", &fDataFilnam);
 263   fPsi2sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 264   fPsi2sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 265   if( fType == 0 ) {
 266     fPsi2sTree ->Branch("fPsi2sESDTracks", &fPsi2sESDTracks);
 267   }
 268   if( fType == 1 ) {
 269     fPsi2sTree ->Branch("fPsi2sAODTracks", &fPsi2sAODTracks);
 270   }
 271   if(isMC) {
 272     fPsi2sTree ->Branch("fGenPart", &fGenPart);
 273     fPsi2sTree ->Branch("fTriggerInputsMC", &fTriggerInputsMC[0], "fTriggerInputsMC[4]/O");
 274   }
 275
 276
 277   fListTrig = new TList();
 278   fListTrig ->SetOwner();
 279
 280   fHistCcup4TriggersPerRun = new TH1D("fHistCcup4TriggersPerRun", "fHistCcup4TriggersPerRun", 33000, 167000.5, 200000.5);
 281   fListTrig->Add(fHistCcup4TriggersPerRun);
 282
 283   fHistCcup7TriggersPerRun = new TH1D("fHistCcup7TriggersPerRun", "fHistCcup7TriggersPerRun", 33000, 167000.5, 200000.5);
 284   fListTrig->Add(fHistCcup7TriggersPerRun);
 285
 286   fHistCcup2TriggersPerRun = new TH1D("fHistCcup2TriggersPerRun", "fHistCcup2TriggersPerRun", 33000, 167000.5, 200000.5);
 287   fListTrig->Add(fHistCcup2TriggersPerRun);
 288
 289   fHistZedTriggersPerRun = new TH1D("fHistZedTriggersPerRun", "fHistZedTriggersPerRun", 33000, 167000.5, 200000.5);
 290   fListTrig->Add(fHistZedTriggersPerRun);
 291
 292   fHistCvlnTriggersPerRun = new TH1D("fHistCvlnTriggersPerRun", "fHistCvlnTriggersPerRun", 33000, 167000.5, 200000.5);
 293   fListTrig->Add(fHistCvlnTriggersPerRun);
 294
 295   fHistMBTriggersPerRun = new TH1D("fHistMBTriggersPerRun", "fHistMBTriggersPerRun", 33000, 167000.5, 200000.5);
 296   fListTrig->Add(fHistMBTriggersPerRun);
 297
 298   fHistCentralTriggersPerRun = new TH1D("fHistCentralTriggersPerRun", "fHistCentralTriggersPerRun", 33000, 167000.5, 200000.5);
 299   fListTrig->Add(fHistCentralTriggersPerRun);
 300
 301   fHistSemiCentralTriggersPerRun = new TH1D("fHistSemiCentralTriggersPerRun", "fHistSemiCentralTriggersPerRun", 33000, 167000.5, 200000.5);
 302   fListTrig->Add(fHistSemiCentralTriggersPerRun);
 303
 304   fListHist = new TList();
 305   fListHist ->SetOwner();
 306
 307   TString CutNameJPsi[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Two good tracks",
 308                                 "Like sign","Oposite sign","One p_{T}>1", "Both p_{T}>1","PID","Dimuom","Dielectron"};
 309   fHistNeventsJPsi = new TH1D("fHistNeventsJPsi","fHistNeventsPsi2s",13,0.5,13.5);
 310   for (Int_t i = 0; i<13; i++) fHistNeventsJPsi->GetXaxis()->SetBinLabel(i+1,CutNameJPsi[i].Data());
 311   fListHist->Add(fHistNeventsJPsi);
 312
 313   fHistTPCsignalJPsi = new TH2D("fHistTPCsignalJPsi","fHistTPCsignalJPsi",240,0,120,240,0,120);
 314   fListHist->Add(fHistTPCsignalJPsi);
 315
 316   fHistDiLeptonPtJPsi = new TH2D("fHistDiLeptonPtJPsi","fHistDiLeptonPtJPsi",350,0,3.5,350,0,3.5);
 317   fListHist->Add(fHistDiLeptonPtJPsi);
 318
 319   fHistDiElectronMass = new TH1D("fHistDiElectronMass","Invariant mass of J/#psi candidates",100,2,5);
 320   fHistDiElectronMass->GetXaxis()->SetTitle("Invariant mass(e^{+}e^{-}) (GeV/c)");
 321   fListHist->Add(fHistDiElectronMass);
 322
 323   fHistDiMuonMass = new TH1D("fHistDiMuonMass","Invariant mass of J/#psi candidates",100,2,5);
 324   fHistDiMuonMass->GetXaxis()->SetTitle("Invariant mass(#mu^{+}#mu^{-}) (GeV/c)");
 325   fListHist->Add(fHistDiMuonMass);
 326
 327   TString CutNamePsi2s[14] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","Four good tracks",
 328                                 "DiLepton - DiPion","Like sign leptons","Like sign pions","Like sign both","Oposite sign","PID","Dimuom","Dielectron"};
 329
 330   fHistNeventsPsi2s = new TH1D("fHistNeventsPsi2s","fHistNeventsPsi2s",14,0.5,14.5);
 331   for (Int_t i = 0; i<14; i++) fHistNeventsPsi2s->GetXaxis()->SetBinLabel(i+1,CutNamePsi2s[i].Data());
 332   fListHist->Add(fHistNeventsPsi2s);
 333
 334   fHistPsi2sMassVsPt = new TH2D("fHistPsi2sMassVsPt","Mass vs p_{T} of #psi(2s) candidates",100,3,6,50,0,5);
 335   fHistPsi2sMassVsPt->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
 336   fHistPsi2sMassVsPt->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 337   fListHist->Add(fHistPsi2sMassVsPt);
 338
 339   fHistPsi2sMassCoherent = new TH1D("fHistPsi2sMassAllCoherent","Invariant mass of coherent #psi(2s) candidates",100,3,6);
 340   fHistPsi2sMassCoherent->GetXaxis()->SetTitle("Invariant mass(l^{+}l^{-}#pi^{+}#pi^{-}) (GeV/c)");
 341   fListHist->Add(fHistPsi2sMassCoherent);
 342
 343   PostData(1, fJPsiTree);
 344   PostData(2, fPsi2sTree);
 345   PostData(3, fListTrig);
 346   PostData(4, fListHist);
 347
 348 }//UserCreateOutputObjects
 349
 350 //_____________________________________________________________________________
 351 void AliAnalysisTaskUpcPsi2s::UserExec(Option_t *)
 352 {
 353
 354   //cout<<"#################### Next event ##################"<<endl;
 355
 356   if( fType == 0 ){
 357         RunESDtrig();
 358         if(fRunHist) RunESDhist();
 359         if(fRunTree) RunESDtree();
 360         }
 361
 362   if( fType == 1 ){
 363         RunAODtrig();
 364         if(fRunHist) RunAODhist();
 365         if(fRunTree) RunAODtree();
 366         }
 367
 368 }//UserExec
 369 //_____________________________________________________________________________
 370 void AliAnalysisTaskUpcPsi2s::RunAODtrig()
 371 {
 372
 373   //input event
 374   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 375   if(!aod) return;
 376
 377   fRunNum = aod ->GetRunNumber();
 378   //Trigger
 379   TString trigger = aod->GetFiredTriggerClasses();
 380
 381   if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
 382   if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
 383   if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
 384
 385   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
 386   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
 387
 388   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 389   if(fL1inputs & (1 << 18)) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
 390
 391   //MB, Central and SemiCentral triggers
 392   AliCentrality *centrality = aod->GetCentrality();
 393   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 394
 395   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
 396   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
 397
 398   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
 399
 400   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
 401
 402   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
 403
 404 PostData(3, fListTrig);
 405
 406 }
 407 //_____________________________________________________________________________
 408 void AliAnalysisTaskUpcPsi2s::RunAODhist()
 409 {
 410
 411   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 412
 413   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 414   Double_t muonMass = partMuon->Mass();
 415
 416   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 417   Double_t electronMass = partElectron->Mass();
 418
 419   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 420   Double_t pionMass = partPion->Mass();
 421
 422   //input event
 423   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 424   if(!aod) return;
 425
 426   fHistNeventsJPsi->Fill(1);
 427   fHistNeventsPsi2s->Fill(1);
 428
 429   //Trigger
 430   TString trigger = aod->GetFiredTriggerClasses();
 431
 432   if(!isMC && !trigger.Contains("CCUP") ) return;
 433
 434   fHistNeventsJPsi->Fill(2);
 435   fHistNeventsPsi2s->Fill(2);
 436
 437   //primary vertex
 438   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 439   fVtxContrib = fAODVertex->GetNContributors();
 440   if(fVtxContrib < 2) return;
 441
 442   fHistNeventsJPsi->Fill(3);
 443   fHistNeventsPsi2s->Fill(3);
 444
 445   //VZERO, ZDC
 446   AliAODVZERO *fV0data = aod ->GetVZEROData();
 447   AliAODZDC *fZDCdata = aod->GetZDCData();
 448
 449   fV0Adecision = fV0data->GetV0ADecision();
 450   fV0Cdecision = fV0data->GetV0CDecision();
 451   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
 452
 453   fHistNeventsJPsi->Fill(4);
 454   fHistNeventsPsi2s->Fill(4);
 455
 456   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 457   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 458
 459   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
 460
 461   fHistNeventsJPsi->Fill(5);
 462   fHistNeventsPsi2s->Fill(5);
 463
 464   //Two tracks loop
 465   Int_t nGoodTracks = 0;
 466   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 467
 468   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
 469   Short_t qLepton[4], qPion[4];
 470   UInt_t nLepton=0, nPion=0, nHighPt=0;
 471   Double_t fRecTPCsignal[5];
 472   Int_t mass[3]={-1,-1,-1};
 473
 474
 475   //Four track loop
 476   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 477     AliAODTrack *trk = aod->GetTrack(itr);
 478     if( !trk ) continue;
 479     if(!(trk->TestFilterBit(1<<0))) continue;
 480
 481       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 482       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 483       if(trk->GetTPCNcls() < 50)continue;
 484       if(trk->Chi2perNDF() > 4)continue;
 485       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 486       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 487       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 488       delete trk_clone;
 489
 490       if(TMath::Abs(dca[1]) > 2) continue;
 491
 492       TrackIndex[nGoodTracks] = itr;
 493       nGoodTracks++;
 494
 495       if(nGoodTracks > 4) break;
 496   }//Track loop
 497
 498   nLepton=0; nPion=0; nHighPt=0;
 499   mass[0]= -1; mass[1]= -1, mass[2]= -1;
 500
 501   if(nGoodTracks == 4){
 502           fHistNeventsPsi2s->Fill(6);
 503           for(Int_t i=0; i<4; i++){
 504                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 505
 506                 if(trk->Pt() > 1){
 507                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 508                         qLepton[nLepton] = trk->Charge();
 509                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 510                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 511                                         mass[nLepton] = 0;
 512                                         }
 513                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 514                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 515                                         mass[nLepton] = 1;
 516                                         }
 517                         nLepton++;
 518                         }
 519                 else{
 520                         qPion[nPion] = trk->Charge();
 521                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
 522                         nPion++;
 523                         }
 524
 525                 if(nLepton > 2 || nPion > 2) break;
 526                 }
 527         if((nLepton == 2) && (nPion == 2)){
 528                 fHistNeventsPsi2s->Fill(7);
 529                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(8);
 530                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(9);
 531                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
 532                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
 533                         fHistNeventsPsi2s->Fill(11);
 534                         if(mass[0] == mass[1]) {
 535                                 fHistNeventsPsi2s->Fill(12);
 536                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
 537                                 vDilepton = vLepton[0]+vLepton[1];
 538                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
 539                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
 540                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(13);
 541                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(14);
 542                                 }
 543                         }
 544                 }
 545   }
 546
 547   nGoodTracks = 0;
 548   //Two track loop
 549   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 550     AliAODTrack *trk = aod->GetTrack(itr);
 551     if( !trk ) continue;
 552     if(!(trk->TestFilterBit(1<<0))) continue;
 553
 554       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 555       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 556       if(trk->GetTPCNcls() < 70)continue;
 557       if(trk->Chi2perNDF() > 4)continue;
 558       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 559       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 560       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 561       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 562       delete trk_clone;
 563       if(TMath::Abs(dca[1]) > 2) continue;
 564       if(TMath::Abs(dca[0]) > 0.2) continue;
 565
 566       TrackIndex[nGoodTracks] = itr;
 567       nGoodTracks++;
 568
 569       if(nGoodTracks > 2) break;
 570   }//Track loop
 571
 572    nLepton=0; nPion=0; nHighPt=0;
 573   mass[0]= -1; mass[1]= -1, mass[2]= -1;
 574
 575   if(nGoodTracks == 2){
 576           fHistNeventsJPsi->Fill(6);
 577           for(Int_t i=0; i<2; i++){
 578                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 579                 if(trk->Pt() > 1) nHighPt++;
 580                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 581                 qLepton[nLepton] = trk->Charge();
 582                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 583                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 584                                 mass[nLepton] = 0;
 585                                 }
 586                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 587                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 588                                 mass[nLepton] = 1;
 589                                 }
 590                 nLepton++;
 591                 }
 592         if(nLepton == 2){
 593                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(7);
 594                 if(qLepton[0]*qLepton[1] < 0){
 595                         fHistNeventsJPsi->Fill(8);
 596                         if(nHighPt > 0){
 597                                 fHistNeventsJPsi->Fill(9);
 598                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
 599                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(10);
 600                                 if(mass[0] == mass[1] && mass[0] != -1) {
 601                                         fHistNeventsJPsi->Fill(11);
 602                                         vCandidate = vLepton[0]+vLepton[1];
 603                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
 604                                         if(mass[0] == 0) {
 605                                                 fHistDiMuonMass->Fill(vCandidate.M());
 606                                                 fHistNeventsJPsi->Fill(12);
 607                                                 }
 608                                         if(mass[0] == 1) {
 609                                                 fHistDiElectronMass->Fill(vCandidate.M());
 610                                                 fHistNeventsJPsi->Fill(13);
 611                                                 }
 612                                         }
 613                                 }
 614                         }
 615                 }
 616   }
 617
 618
 619   PostData(4, fListHist);
 620
 621 }
 622
 623 //_____________________________________________________________________________
 624 void AliAnalysisTaskUpcPsi2s::RunAODtree()
 625 {
 626   //input event
 627   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 628   if(!aod) return;
 629
 630   if(isMC) RunAODMC(aod);
 631
 632   //input data
 633   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
 634   fDataFilnam->Clear();
 635   fDataFilnam->SetString(filnam);
 636   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
 637   fRunNum = aod ->GetRunNumber();
 638
 639   //Trigger
 640   TString trigger = aod->GetFiredTriggerClasses();
 641
 642   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
 643   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
 644   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
 645
 646   Bool_t isTriggered = kFALSE;
 647   for(Int_t i=0; i<ntrg; i++) {
 648     if( fTrigger[i] ) isTriggered = kTRUE;
 649   }
 650   if(!isMC && !isTriggered ) return;
 651
 652   //trigger inputs
 653   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
 654   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 655
 656   //Event identification
 657   fPerNum = aod ->GetPeriodNumber();
 658   fOrbNum = aod ->GetOrbitNumber();
 659   fBCrossNum = aod ->GetBunchCrossNumber();
 660
 661   //primary vertex
 662   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 663   fVtxContrib = fAODVertex->GetNContributors();
 664   fVtxPos[0] = fAODVertex->GetX();
 665   fVtxPos[1] = fAODVertex->GetY();
 666   fVtxPos[2] = fAODVertex->GetZ();
 667   Double_t CovMatx[6];
 668   fAODVertex->GetCovarianceMatrix(CovMatx);
 669   fVtxErr[0] = CovMatx[0];
 670   fVtxErr[1] = CovMatx[1];
 671   fVtxErr[2] = CovMatx[2];
 672   fVtxChi2 = fAODVertex->GetChi2();
 673   fVtxNDF = fAODVertex->GetNDF();
 674
 675   //Tracklets
 676   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
 677
 678   //VZERO, ZDC
 679   AliAODVZERO *fV0data = aod ->GetVZEROData();
 680   AliAODZDC *fZDCdata = aod->GetZDCData();
 681
 682   fV0Adecision = fV0data->GetV0ADecision();
 683   fV0Cdecision = fV0data->GetV0CDecision();
 684   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 685   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 686
 687   fNLooseTracks = 0;
 688
 689   //Track loop - loose cuts
 690   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 691     AliAODTrack *trk = aod->GetTrack(itr);
 692     if( !trk ) continue;
 693     if(!(trk->TestFilterBit(1<<0))) continue;
 694
 695       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 696       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 697       if(trk->GetTPCNcls() < 20)continue;
 698       fNLooseTracks++;
 699   }//Track loop -loose cuts
 700
 701   Int_t nGoodTracks=0;
 702   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 703
 704   //Two track loop
 705   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 706     AliAODTrack *trk = aod->GetTrack(itr);
 707     if( !trk ) continue;
 708     if(!(trk->TestFilterBit(1<<0))) continue;
 709
 710       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 711       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 712       if(trk->GetTPCNcls() < 70)continue;
 713       if(trk->Chi2perNDF() > 4)continue;
 714       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 715       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 716       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 717       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 718       delete trk_clone;
 719       if(TMath::Abs(dca[1]) > 2) continue;
 720       Double_t cut_DCAxy = (0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 721       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
 722
 723
 724       TrackIndex[nGoodTracks] = itr;
 725       nGoodTracks++;
 726
 727       if(nGoodTracks > 2) break;
 728   }//Track loop
 729
 730   fJPsiAODTracks->Clear("C");
 731   if(nGoodTracks == 2){
 732
 733           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 734           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 735           Double_t muonMass = partMuon->Mass();
 736           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 737           Double_t electronMass = partElectron->Mass();
 738           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 739           Double_t pionMass = partPion->Mass();
 740
 741           Double_t KFcov[21];
 742           Double_t KFpar[6];
 743           Double_t KFmass = pionMass;
 744           Double_t fRecTPCsignal;
 745           AliKFParticle *KFpart[2];
 746           AliKFVertex *KFvtx = new AliKFVertex();
 747           KFvtx->SetField(aod->GetMagneticField());
 748
 749           for(Int_t i=0; i<2; i++){
 750                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 751
 752                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 753                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 754                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 755                 delete trk_clone;
 756
 757                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
 758                 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
 759
 760                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
 761                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
 762                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
 763
 764                 trk->GetPosition(KFpar);
 765                 trk->PxPyPz(KFpar+3);
 766                 trk->GetCovMatrix(KFcov);
 767
 768                 if(trk->Pt() > 1){
 769                         fRecTPCsignal = trk->GetTPCsignal();
 770                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
 771                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
 772                         }
 773                 else KFmass = pionMass;
 774
 775                 KFpart[i] = new AliKFParticle();
 776                 KFpart[i]->SetField(aod->GetMagneticField());
 777                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
 778                 KFvtx->AddDaughter(*KFpart[i]);
 779
 780
 781                 Double_t pos[3]={0,0,0};
 782                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
 783                 parTrk->CopyFromVTrack((AliVTrack*) trk);
 784                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
 785                 else {
 786                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
 787                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
 788                      }
 789                 delete parTrk;
 790                 }
 791   fKfVtxPos[0]= KFvtx->GetX();
 792   fKfVtxPos[1]= KFvtx->GetY();
 793   fKfVtxPos[2]= KFvtx->GetZ();
 794   for(UInt_t i=0; i<2; i++)delete KFpart[i];
 795   delete KFvtx;
 796
 797   if(!isMC) fJPsiTree ->Fill();
 798   }
 799
 800    nGoodTracks = 0;
 801    //Four track loop
 802   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 803     AliAODTrack *trk = aod->GetTrack(itr);
 804     if( !trk ) continue;
 805     if(!(trk->TestFilterBit(1<<0))) continue;
 806
 807       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 808       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 809       if(trk->GetTPCNcls() < 50)continue;
 810       if(trk->Chi2perNDF() > 4)continue;
 811       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 812       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 813       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 814       delete trk_clone;
 815       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 816       if(TMath::Abs(dca[1]) > 2) continue;
 817       Double_t cut_DCAxy = 4*(0.0182 + 0.0350/TMath::Power(trk->Pt(),1.01));
 818       if(TMath::Abs(dca[0]) > cut_DCAxy) continue;
 819
 820       TrackIndex[nGoodTracks] = itr;
 821       nGoodTracks++;
 822
 823       if(nGoodTracks > 4) break;
 824   }//Track loop
 825
 826   fPsi2sAODTracks->Clear("C");
 827   if(nGoodTracks == 4){
 828
 829           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 830           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 831           Double_t muonMass = partMuon->Mass();
 832           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 833           Double_t electronMass = partElectron->Mass();
 834           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 835           Double_t pionMass = partPion->Mass();
 836
 837           Double_t KFcov[21];
 838           Double_t KFpar[6];
 839           Double_t KFmass = pionMass;
 840           Double_t fRecTPCsignal;
 841           AliKFParticle *KFpart[4];
 842           AliKFVertex *KFvtx = new AliKFVertex();
 843           KFvtx->SetField(aod->GetMagneticField());
 844
 845           for(Int_t i=0; i<4; i++){
 846                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 847
 848                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 849                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 850                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 851                 delete trk_clone;
 852
 853                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
 854                 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
 855
 856
 857                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
 858                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
 859                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
 860
 861                 trk->GetPosition(KFpar);
 862                 trk->PxPyPz(KFpar+3);
 863                 trk->GetCovMatrix(KFcov);
 864
 865                 if(trk->Pt() > 1){
 866                         fRecTPCsignal = trk->GetTPCsignal();
 867                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
 868                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
 869                         }
 870                 else KFmass = pionMass;
 871
 872                 KFpart[i] = new AliKFParticle();
 873                 KFpart[i]->SetField(aod->GetMagneticField());
 874                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
 875                 KFvtx->AddDaughter(*KFpart[i]);
 876
 877                 Double_t pos[3]={0,0,0};
 878                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
 879                 parTrk->CopyFromVTrack((AliVTrack*) trk);
 880                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
 881                 else {
 882                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
 883                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
 884                      }
 885                 delete parTrk;
 886                 }
 887   fKfVtxPos[0]= KFvtx->GetX();
 888   fKfVtxPos[1]= KFvtx->GetY();
 889   fKfVtxPos[2]= KFvtx->GetZ();
 890   for(UInt_t i=0; i<4; i++)delete KFpart[i];
 891   delete KFvtx;
 892   if(!isMC) fPsi2sTree ->Fill();
 893   }
 894
 895   if(isMC){
 896         fJPsiTree ->Fill();
 897         fPsi2sTree ->Fill();
 898   }
 899
 900   PostData(1, fJPsiTree);
 901   PostData(2, fPsi2sTree);
 902
 903 }//RunAOD
 904
 905
 906 //_____________________________________________________________________________
 907 void AliAnalysisTaskUpcPsi2s::RunAODMC(AliAODEvent *aod)
 908 {
 909
 910   fGenPart->Clear("C");
 911
 912   TClonesArray *arrayMC = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 913   if(!arrayMC) return;
 914
 915   Int_t nmc=0;
 916   //loop over mc particles
 917   for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
 918     AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(imc);
 919     if(!mcPart) continue;
 920
 921     if(mcPart->GetMother() >= 0) continue;
 922
 923     TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
 924     part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
 925     part->SetPdgCode(mcPart->GetPdgCode());
 926     part->SetUniqueID(imc);
 927   }//loop over mc particles
 928
 929 }//RunAODMC
 930
 931
 932 //_____________________________________________________________________________
 933 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
 934 {
 935
 936   //input event
 937   AliESDEvent *esd = (AliESDEvent*) InputEvent();
 938   if(!esd) return;
 939
 940   fRunNum = esd ->GetRunNumber();
 941   //Trigger
 942   TString trigger = esd->GetFiredTriggerClasses();
 943
 944   if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
 945   if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
 946   if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
 947
 948   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
 949   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
 950
 951   if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
 952
 953    //MB, Central and SemiCentral triggers
 954   AliCentrality *centrality = esd->GetCentrality();
 955   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 956
 957   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
 958   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
 959
 960   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
 961
 962   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
 963
 964   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
 965
 966
 967 PostData(3, fListTrig);
 968
 969 }
 970 //_____________________________________________________________________________
 971 void AliAnalysisTaskUpcPsi2s::RunESDhist()
 972 {
 973
 974   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 975
 976   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 977   Double_t muonMass = partMuon->Mass();
 978
 979   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 980   Double_t electronMass = partElectron->Mass();
 981
 982   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 983   Double_t pionMass = partPion->Mass();
 984
 985   //input event
 986   AliESDEvent *esd = (AliESDEvent*) InputEvent();
 987   if(!esd) return;
 988
 989   fHistNeventsJPsi->Fill(1);
 990   fHistNeventsPsi2s->Fill(1);
 991
 992   //Trigger
 993   TString trigger = esd->GetFiredTriggerClasses();
 994
 995   if(!isMC && !trigger.Contains("CCUP") ) return;
 996
 997   fHistNeventsJPsi->Fill(2);
 998   fHistNeventsPsi2s->Fill(2);
 999
1000   //primary vertex
1001   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1002   fVtxContrib = fESDVertex->GetNContributors();
1003   if(fVtxContrib < 2) return;
1004
1005   fHistNeventsJPsi->Fill(3);
1006   fHistNeventsPsi2s->Fill(3);
1007
1008   //VZERO, ZDC
1009   AliESDVZERO *fV0data = esd->GetVZEROData();
1010   AliESDZDC *fZDCdata = esd->GetESDZDC();
1011
1012   fV0Adecision = fV0data->GetV0ADecision();
1013   fV0Cdecision = fV0data->GetV0CDecision();
1014   if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
1015
1016   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1017   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1018   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
1019
1020   fHistNeventsJPsi->Fill(4);
1021   fHistNeventsPsi2s->Fill(4);
1022
1023    Int_t nGoodTracks=0;
1024   //Two tracks loop
1025   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1026
1027   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1028   Short_t qLepton[4], qPion[4];
1029   UInt_t nLepton=0, nPion=0, nHighPt=0;
1030   Double_t fRecTPCsignal[5];
1031   Int_t mass[3]={-1,-1,-1};
1032
1033  //Two Track loop
1034   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1035     AliESDtrack *trk = esd->GetTrack(itr);
1036     if( !trk ) continue;
1037
1038       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1039       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1040       if(trk->GetTPCNcls() < 70)continue;
1041       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1042       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1043       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1044       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1045       trk->GetImpactParameters(dca[0],dca[1]);
1046       if(TMath::Abs(dca[1]) > 2) continue;
1047       if(TMath::Abs(dca[1]) > 0.2) continue;
1048
1049       TrackIndex[nGoodTracks] = itr;
1050       nGoodTracks++;
1051       if(nGoodTracks > 2) break;
1052   }//Track loop
1053
1054
1055   if(nGoodTracks == 2){
1056           fHistNeventsJPsi->Fill(5);
1057           for(Int_t i=0; i<2; i++){
1058                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1059                 if(trk->Pt() > 1) nHighPt++;
1060                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1061                 qLepton[nLepton] = trk->Charge();
1062                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1063                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1064                                 mass[nLepton] = 0;
1065                                 }
1066                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1067                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1068                                 mass[nLepton] = 1;
1069                                 }
1070                 nLepton++;
1071                 }
1072         if(nLepton == 2){
1073                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
1074                 if(qLepton[0]*qLepton[1] < 0){
1075                         fHistNeventsJPsi->Fill(7);
1076                         if(nHighPt > 0){
1077                                 fHistNeventsJPsi->Fill(8);
1078                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
1079                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
1080                                 if(mass[0] == mass[1] && mass[0] != -1) {
1081                                         fHistNeventsJPsi->Fill(10);
1082                                         vCandidate = vLepton[0]+vLepton[1];
1083                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
1084                                         if(mass[0] == 0) {
1085                                                 fHistDiMuonMass->Fill(vCandidate.M());
1086                                                 fHistNeventsJPsi->Fill(11);
1087                                                 }
1088                                         if(mass[0] == 1) {
1089                                                 fHistDiElectronMass->Fill(vCandidate.M());
1090                                                 fHistNeventsJPsi->Fill(12);
1091                                                 }
1092                                         }
1093                                 }
1094                         }
1095                 }
1096   }
1097   nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1098   mass[0]= -1; mass[1]= -1, mass[2]= -1;
1099
1100     //Four Track loop
1101   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1102     AliESDtrack *trk = esd->GetTrack(itr);
1103     if( !trk ) continue;
1104
1105       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1106       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1107       if(trk->GetTPCNcls() < 50)continue;
1108       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1109       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1110       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1111       trk->GetImpactParameters(dca[0],dca[1]);
1112       if(TMath::Abs(dca[1]) > 2) continue;
1113
1114       TrackIndex[nGoodTracks] = itr;
1115       nGoodTracks++;
1116       if(nGoodTracks > 4) break;
1117   }//Track loop
1118
1119   if(nGoodTracks == 4){
1120           fHistNeventsPsi2s->Fill(5);
1121           for(Int_t i=0; i<4; i++){
1122                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1123
1124                 if(trk->Pt() > 1){
1125                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1126                         qLepton[nLepton] = trk->Charge();
1127                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1128                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1129                                         mass[nLepton] = 0;
1130                                         }
1131                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1132                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1133                                         mass[nLepton] = 1;
1134                                         }
1135                         nLepton++;
1136                         }
1137                 else{
1138                         qPion[nPion] = trk->Charge();
1139                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1140                         nPion++;
1141                         }
1142
1143                 if(nLepton > 2 || nPion > 2) break;
1144                 }
1145         if((nLepton == 2) && (nPion == 2)){
1146                 fHistNeventsPsi2s->Fill(6);
1147                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1148                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1149                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1150                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1151                         fHistNeventsPsi2s->Fill(10);
1152                         if(mass[0] == mass[1]) {
1153                                 fHistNeventsPsi2s->Fill(11);
1154                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1155                                 vDilepton = vLepton[0]+vLepton[1];
1156                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1157                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1158                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1159                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1160                                 }
1161                         }
1162                 }
1163   }
1164
1165   PostData(4, fListHist);
1166
1167 }
1168
1169 //_____________________________________________________________________________
1170 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1171 {
1172
1173   //input event
1174   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1175   if(!esd) return;
1176
1177   if(isMC) RunESDMC(esd);
1178
1179   //input data
1180   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1181   fDataFilnam->Clear();
1182   fDataFilnam->SetString(filnam);
1183   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1184   fRunNum = esd->GetRunNumber();
1185
1186    //Trigger
1187   TString trigger = esd->GetFiredTriggerClasses();
1188
1189   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1190   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
1191   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1192
1193   Bool_t isTriggered = kFALSE;
1194   for(Int_t i=0; i<ntrg; i++) {
1195     if( fTrigger[i] ) isTriggered = kTRUE;
1196   }
1197   if(!isMC && !isTriggered ) return;
1198
1199   //trigger inputs
1200   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1201   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1202
1203   //TOF trigger info (0OMU)
1204   fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1205   fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1206
1207   //Event identification
1208   fPerNum = esd->GetPeriodNumber();
1209   fOrbNum = esd->GetOrbitNumber();
1210   fBCrossNum = esd->GetBunchCrossNumber();
1211
1212   //primary vertex
1213   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1214   fVtxContrib = fESDVertex->GetNContributors();
1215   fVtxPos[0] = fESDVertex->GetX();
1216   fVtxPos[1] = fESDVertex->GetY();
1217   fVtxPos[2] = fESDVertex->GetZ();
1218   Double_t CovMatx[6];
1219   fESDVertex->GetCovarianceMatrix(CovMatx);
1220   fVtxErr[0] = CovMatx[0];
1221   fVtxErr[1] = CovMatx[1];
1222   fVtxErr[2] = CovMatx[2];
1223   fVtxChi2 = fESDVertex->GetChi2();
1224   fVtxNDF = fESDVertex->GetNDF();
1225
1226   //Tracklets
1227   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1228
1229   //VZERO, ZDC
1230   AliESDVZERO *fV0data = esd->GetVZEROData();
1231   AliESDZDC *fZDCdata = esd->GetESDZDC();
1232
1233   fV0Adecision = fV0data->GetV0ADecision();
1234   fV0Cdecision = fV0data->GetV0CDecision();
1235   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1236   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1237
1238   fNLooseTracks = 0;
1239
1240   //Track loop - loose cuts
1241   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1242     AliESDtrack *trk = esd->GetTrack(itr);
1243     if( !trk ) continue;
1244
1245       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1246       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1247       if(trk->GetTPCNcls() < 20)continue;
1248       fNLooseTracks++;
1249   }//Track loop -loose cuts
1250
1251   Int_t nGoodTracks=0;
1252   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1253
1254   //Two Track loop
1255   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1256     AliESDtrack *trk = esd->GetTrack(itr);
1257     if( !trk ) continue;
1258
1259       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1260       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1261       if(trk->GetTPCNcls() < 70)continue;
1262       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1263       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1264       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1265       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1266       trk->GetImpactParameters(dca[0],dca[1]);
1267       if(TMath::Abs(dca[1]) > 2) continue;
1268       if(TMath::Abs(dca[1]) > 0.2) continue;
1269
1270       TrackIndex[nGoodTracks] = itr;
1271       nGoodTracks++;
1272       if(nGoodTracks > 2) break;
1273   }//Track loop
1274
1275   fJPsiESDTracks->Clear("C");
1276   if(nGoodTracks == 2){
1277           for(Int_t i=0; i<2; i++){
1278                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1279
1280                 AliExternalTrackParam cParam;
1281                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1282
1283                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
1284
1285                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1286                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1287                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1288
1289                 Double_t pos[3]={0,0,0};
1290                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1291                 else {
1292                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1293                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1294                      }
1295                 }
1296   if(!isMC) fJPsiTree ->Fill();
1297   }
1298
1299   nGoodTracks = 0;
1300   //Four track loop
1301   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1302     AliESDtrack *trk = esd->GetTrack(itr);
1303     if( !trk ) continue;
1304
1305       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1306       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1307       if(trk->GetTPCNcls() < 50)continue;
1308       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1309       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1310       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1311       trk->GetImpactParameters(dca[0],dca[1]);
1312       if(TMath::Abs(dca[1]) > 2) continue;
1313
1314       TrackIndex[nGoodTracks] = itr;
1315       nGoodTracks++;
1316       if(nGoodTracks > 4) break;
1317   }//Track loop
1318
1319   fPsi2sESDTracks->Clear("C");
1320   if(nGoodTracks == 4){
1321           for(Int_t i=0; i<4; i++){
1322                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1323
1324                 AliExternalTrackParam cParam;
1325                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1326
1327                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1328
1329                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1330                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1331                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1332
1333                 Double_t pos[3]={0,0,0};
1334                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1335                 else {
1336                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1337                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1338                      }
1339                 }
1340   if(!isMC) fPsi2sTree ->Fill();
1341   }
1342
1343   if(isMC){
1344         fJPsiTree ->Fill();
1345         fPsi2sTree ->Fill();
1346   }
1347
1348   PostData(1, fJPsiTree);
1349   PostData(2, fPsi2sTree);
1350
1351 }//RunESD
1352
1353
1354 //_____________________________________________________________________________
1355 void AliAnalysisTaskUpcPsi2s::RunESDMC(AliESDEvent* esd)
1356 {
1357
1358   AliTriggerAnalysis *fTrigAna = new AliTriggerAnalysis();
1359   fTrigAna->SetAnalyzeMC(isMC);
1360
1361   if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) > 1) fTriggerInputsMC[0] = kTRUE;
1362   if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) < 2) fTriggerInputsMC[0] = kFALSE;
1363
1364   fTriggerInputsMC[1] = esd->GetHeader()->IsTriggerInputFired("0OMU");
1365   fTriggerInputsMC[2] = esd->GetHeader()->IsTriggerInputFired("0VBA");
1366   fTriggerInputsMC[3] = esd->GetHeader()->IsTriggerInputFired("0VBC");
1367
1368   fGenPart->Clear("C");
1369
1370   AliMCEvent *mc = MCEvent();
1371   if(!mc) return;
1372
1373   Int_t nmc = 0;
1374   //loop over mc particles
1375   for(Int_t imc=0; imc<mc->GetNumberOfTracks(); imc++) {
1376     AliMCParticle *mcPart = (AliMCParticle*) mc->GetTrack(imc);
1377     if(!mcPart) continue;
1378
1379     if(mcPart->GetMother() >= 0) continue;
1380
1381     TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1382     part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1383     part->SetPdgCode(mcPart->PdgCode());
1384     part->SetUniqueID(imc);
1385   }//loop over mc particles
1386
1387 }//RunESDMC
1388
1389
1390
1391 //_____________________________________________________________________________
1392 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
1393 {
1394
1395   cout<<"Analysis complete."<<endl;
1396 }//Terminate
1397