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       if(TMath::Abs(dca[0]) > 0.2) continue;
 721
 722       TrackIndex[nGoodTracks] = itr;
 723       nGoodTracks++;
 724
 725       if(nGoodTracks > 2) break;
 726   }//Track loop
 727
 728   fJPsiAODTracks->Clear("C");
 729   if(nGoodTracks == 2){
 730
 731           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 732           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 733           Double_t muonMass = partMuon->Mass();
 734           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 735           Double_t electronMass = partElectron->Mass();
 736           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 737           Double_t pionMass = partPion->Mass();
 738
 739           Double_t KFcov[21];
 740           Double_t KFpar[6];
 741           Double_t KFmass = pionMass;
 742           Double_t fRecTPCsignal;
 743           AliKFParticle *KFpart[2];
 744           AliKFVertex *KFvtx = new AliKFVertex();
 745           KFvtx->SetField(aod->GetMagneticField());
 746
 747           for(Int_t i=0; i<2; i++){
 748                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 749
 750                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 751                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 752                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 753                 delete trk_clone;
 754
 755                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
 756                 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
 757
 758                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
 759                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
 760                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
 761
 762                 trk->GetPosition(KFpar);
 763                 trk->PxPyPz(KFpar+3);
 764                 trk->GetCovMatrix(KFcov);
 765
 766                 if(trk->Pt() > 1){
 767                         fRecTPCsignal = trk->GetTPCsignal();
 768                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
 769                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
 770                         }
 771                 else KFmass = pionMass;
 772
 773                 KFpart[i] = new AliKFParticle();
 774                 KFpart[i]->SetField(aod->GetMagneticField());
 775                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
 776                 KFvtx->AddDaughter(*KFpart[i]);
 777
 778
 779                 Double_t pos[3]={0,0,0};
 780                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
 781                 parTrk->CopyFromVTrack((AliVTrack*) trk);
 782                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
 783                 else {
 784                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
 785                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
 786                      }
 787                 delete parTrk;
 788                 }
 789   fKfVtxPos[0]= KFvtx->GetX();
 790   fKfVtxPos[1]= KFvtx->GetY();
 791   fKfVtxPos[2]= KFvtx->GetZ();
 792   for(UInt_t i=0; i<2; i++)delete KFpart[i];
 793   delete KFvtx;
 794
 795   if(!isMC) fJPsiTree ->Fill();
 796   }
 797
 798    nGoodTracks = 0;
 799    //Four track loop
 800   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 801     AliAODTrack *trk = aod->GetTrack(itr);
 802     if( !trk ) continue;
 803     if(!(trk->TestFilterBit(1<<0))) continue;
 804
 805       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 806       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 807       if(trk->GetTPCNcls() < 50)continue;
 808       if(trk->Chi2perNDF() > 4)continue;
 809       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 810       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 811       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 812       delete trk_clone;
 813       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 814       if(TMath::Abs(dca[1]) > 2) continue;
 815
 816       TrackIndex[nGoodTracks] = itr;
 817       nGoodTracks++;
 818
 819       if(nGoodTracks > 4) break;
 820   }//Track loop
 821
 822   fPsi2sAODTracks->Clear("C");
 823   if(nGoodTracks == 4){
 824
 825           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 826           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 827           Double_t muonMass = partMuon->Mass();
 828           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 829           Double_t electronMass = partElectron->Mass();
 830           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 831           Double_t pionMass = partPion->Mass();
 832
 833           Double_t KFcov[21];
 834           Double_t KFpar[6];
 835           Double_t KFmass = pionMass;
 836           Double_t fRecTPCsignal;
 837           AliKFParticle *KFpart[4];
 838           AliKFVertex *KFvtx = new AliKFVertex();
 839           KFvtx->SetField(aod->GetMagneticField());
 840
 841           for(Int_t i=0; i<4; i++){
 842                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 843
 844                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 845                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 846                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 847                 delete trk_clone;
 848
 849                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
 850                 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
 851
 852
 853                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
 854                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
 855                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
 856
 857                 trk->GetPosition(KFpar);
 858                 trk->PxPyPz(KFpar+3);
 859                 trk->GetCovMatrix(KFcov);
 860
 861                 if(trk->Pt() > 1){
 862                         fRecTPCsignal = trk->GetTPCsignal();
 863                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
 864                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
 865                         }
 866                 else KFmass = pionMass;
 867
 868                 KFpart[i] = new AliKFParticle();
 869                 KFpart[i]->SetField(aod->GetMagneticField());
 870                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
 871                 KFvtx->AddDaughter(*KFpart[i]);
 872
 873                 Double_t pos[3]={0,0,0};
 874                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
 875                 parTrk->CopyFromVTrack((AliVTrack*) trk);
 876                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
 877                 else {
 878                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
 879                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
 880                      }
 881                 delete parTrk;
 882                 }
 883   fKfVtxPos[0]= KFvtx->GetX();
 884   fKfVtxPos[1]= KFvtx->GetY();
 885   fKfVtxPos[2]= KFvtx->GetZ();
 886   for(UInt_t i=0; i<4; i++)delete KFpart[i];
 887   delete KFvtx;
 888   if(!isMC) fPsi2sTree ->Fill();
 889   }
 890
 891   if(isMC){
 892         fJPsiTree ->Fill();
 893         fPsi2sTree ->Fill();
 894   }
 895
 896   PostData(1, fJPsiTree);
 897   PostData(2, fPsi2sTree);
 898
 899 }//RunAOD
 900
 901
 902 //_____________________________________________________________________________
 903 void AliAnalysisTaskUpcPsi2s::RunAODMC(AliAODEvent *aod)
 904 {
 905
 906   fGenPart->Clear("C");
 907
 908   TClonesArray *arrayMC = (TClonesArray*) aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 909   if(!arrayMC) return;
 910
 911   Int_t nmc=0;
 912   //loop over mc particles
 913   for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
 914     AliAODMCParticle *mcPart = (AliAODMCParticle*) arrayMC->At(imc);
 915     if(!mcPart) continue;
 916
 917     if(mcPart->GetMother() >= 0) continue;
 918
 919     TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
 920     part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
 921     part->SetPdgCode(mcPart->GetPdgCode());
 922     part->SetUniqueID(imc);
 923   }//loop over mc particles
 924
 925 }//RunAODMC
 926
 927
 928 //_____________________________________________________________________________
 929 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
 930 {
 931
 932   //input event
 933   AliESDEvent *esd = (AliESDEvent*) InputEvent();
 934   if(!esd) return;
 935
 936   fRunNum = esd ->GetRunNumber();
 937   //Trigger
 938   TString trigger = esd->GetFiredTriggerClasses();
 939
 940   if(trigger.Contains("CCUP4-B")) fHistCcup4TriggersPerRun->Fill(fRunNum); //CCUP4 triggers
 941   if(trigger.Contains("CCUP7-B")) fHistCcup7TriggersPerRun->Fill(fRunNum); //CCUP7 triggers
 942   if(trigger.Contains("CCUP2-B")) fHistCcup2TriggersPerRun->Fill(fRunNum); //CCUP2 triggers
 943
 944   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
 945   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
 946
 947   if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
 948
 949    //MB, Central and SemiCentral triggers
 950   AliCentrality *centrality = esd->GetCentrality();
 951   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 952
 953   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
 954   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
 955
 956   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
 957
 958   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0 && (trigger.Contains("CVHN_R2-B"))) fHistCentralTriggersPerRun->Fill(fRunNum);
 959
 960   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
 961
 962
 963 PostData(3, fListTrig);
 964
 965 }
 966 //_____________________________________________________________________________
 967 void AliAnalysisTaskUpcPsi2s::RunESDhist()
 968 {
 969
 970   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 971
 972   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 973   Double_t muonMass = partMuon->Mass();
 974
 975   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 976   Double_t electronMass = partElectron->Mass();
 977
 978   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 979   Double_t pionMass = partPion->Mass();
 980
 981   //input event
 982   AliESDEvent *esd = (AliESDEvent*) InputEvent();
 983   if(!esd) return;
 984
 985   fHistNeventsJPsi->Fill(1);
 986   fHistNeventsPsi2s->Fill(1);
 987
 988   //Trigger
 989   TString trigger = esd->GetFiredTriggerClasses();
 990
 991   if(!isMC && !trigger.Contains("CCUP") ) return;
 992
 993   fHistNeventsJPsi->Fill(2);
 994   fHistNeventsPsi2s->Fill(2);
 995
 996   //primary vertex
 997   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
 998   fVtxContrib = fESDVertex->GetNContributors();
 999   if(fVtxContrib < 2) return;
1000
1001   fHistNeventsJPsi->Fill(3);
1002   fHistNeventsPsi2s->Fill(3);
1003
1004   //VZERO, ZDC
1005   AliESDVZERO *fV0data = esd->GetVZEROData();
1006   AliESDZDC *fZDCdata = esd->GetESDZDC();
1007
1008   fV0Adecision = fV0data->GetV0ADecision();
1009   fV0Cdecision = fV0data->GetV0CDecision();
1010   if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
1011
1012   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1013   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1014   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
1015
1016   fHistNeventsJPsi->Fill(4);
1017   fHistNeventsPsi2s->Fill(4);
1018
1019    Int_t nGoodTracks=0;
1020   //Two tracks loop
1021   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1022
1023   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
1024   Short_t qLepton[4], qPion[4];
1025   UInt_t nLepton=0, nPion=0, nHighPt=0;
1026   Double_t fRecTPCsignal[5];
1027   Int_t mass[3]={-1,-1,-1};
1028
1029  //Two Track loop
1030   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1031     AliESDtrack *trk = esd->GetTrack(itr);
1032     if( !trk ) continue;
1033
1034       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1035       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1036       if(trk->GetTPCNcls() < 70)continue;
1037       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1038       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1039       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1040       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1041       trk->GetImpactParameters(dca[0],dca[1]);
1042       if(TMath::Abs(dca[1]) > 2) continue;
1043       if(TMath::Abs(dca[1]) > 0.2) continue;
1044
1045       TrackIndex[nGoodTracks] = itr;
1046       nGoodTracks++;
1047       if(nGoodTracks > 2) break;
1048   }//Track loop
1049
1050
1051   if(nGoodTracks == 2){
1052           fHistNeventsJPsi->Fill(5);
1053           for(Int_t i=0; i<2; i++){
1054                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1055                 if(trk->Pt() > 1) nHighPt++;
1056                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1057                 qLepton[nLepton] = trk->Charge();
1058                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1059                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1060                                 mass[nLepton] = 0;
1061                                 }
1062                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1063                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1064                                 mass[nLepton] = 1;
1065                                 }
1066                 nLepton++;
1067                 }
1068         if(nLepton == 2){
1069                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
1070                 if(qLepton[0]*qLepton[1] < 0){
1071                         fHistNeventsJPsi->Fill(7);
1072                         if(nHighPt > 0){
1073                                 fHistNeventsJPsi->Fill(8);
1074                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
1075                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
1076                                 if(mass[0] == mass[1] && mass[0] != -1) {
1077                                         fHistNeventsJPsi->Fill(10);
1078                                         vCandidate = vLepton[0]+vLepton[1];
1079                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
1080                                         if(mass[0] == 0) {
1081                                                 fHistDiMuonMass->Fill(vCandidate.M());
1082                                                 fHistNeventsJPsi->Fill(11);
1083                                                 }
1084                                         if(mass[0] == 1) {
1085                                                 fHistDiElectronMass->Fill(vCandidate.M());
1086                                                 fHistNeventsJPsi->Fill(12);
1087                                                 }
1088                                         }
1089                                 }
1090                         }
1091                 }
1092   }
1093   nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1094   mass[0]= -1; mass[1]= -1, mass[2]= -1;
1095
1096     //Four Track loop
1097   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1098     AliESDtrack *trk = esd->GetTrack(itr);
1099     if( !trk ) continue;
1100
1101       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1102       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1103       if(trk->GetTPCNcls() < 50)continue;
1104       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1105       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1106       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1107       trk->GetImpactParameters(dca[0],dca[1]);
1108       if(TMath::Abs(dca[1]) > 2) continue;
1109
1110       TrackIndex[nGoodTracks] = itr;
1111       nGoodTracks++;
1112       if(nGoodTracks > 4) break;
1113   }//Track loop
1114
1115   if(nGoodTracks == 4){
1116           fHistNeventsPsi2s->Fill(5);
1117           for(Int_t i=0; i<4; i++){
1118                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1119
1120                 if(trk->Pt() > 1){
1121                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1122                         qLepton[nLepton] = trk->Charge();
1123                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1124                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1125                                         mass[nLepton] = 0;
1126                                         }
1127                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1128                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1129                                         mass[nLepton] = 1;
1130                                         }
1131                         nLepton++;
1132                         }
1133                 else{
1134                         qPion[nPion] = trk->Charge();
1135                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1136                         nPion++;
1137                         }
1138
1139                 if(nLepton > 2 || nPion > 2) break;
1140                 }
1141         if((nLepton == 2) && (nPion == 2)){
1142                 fHistNeventsPsi2s->Fill(6);
1143                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1144                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1145                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1146                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1147                         fHistNeventsPsi2s->Fill(10);
1148                         if(mass[0] == mass[1]) {
1149                                 fHistNeventsPsi2s->Fill(11);
1150                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1151                                 vDilepton = vLepton[0]+vLepton[1];
1152                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1153                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1154                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1155                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1156                                 }
1157                         }
1158                 }
1159   }
1160
1161   PostData(4, fListHist);
1162
1163 }
1164
1165 //_____________________________________________________________________________
1166 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1167 {
1168
1169   //input event
1170   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1171   if(!esd) return;
1172
1173   if(isMC) RunESDMC(esd);
1174
1175   //input data
1176   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1177   fDataFilnam->Clear();
1178   fDataFilnam->SetString(filnam);
1179   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1180   fRunNum = esd->GetRunNumber();
1181
1182    //Trigger
1183   TString trigger = esd->GetFiredTriggerClasses();
1184
1185   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1186   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
1187   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1188
1189   Bool_t isTriggered = kFALSE;
1190   for(Int_t i=0; i<ntrg; i++) {
1191     if( fTrigger[i] ) isTriggered = kTRUE;
1192   }
1193   if(!isMC && !isTriggered ) return;
1194
1195   //trigger inputs
1196   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1197   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1198
1199   //TOF trigger info (0OMU)
1200   fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1201   fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1202
1203   //Event identification
1204   fPerNum = esd->GetPeriodNumber();
1205   fOrbNum = esd->GetOrbitNumber();
1206   fBCrossNum = esd->GetBunchCrossNumber();
1207
1208   //primary vertex
1209   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1210   fVtxContrib = fESDVertex->GetNContributors();
1211   fVtxPos[0] = fESDVertex->GetX();
1212   fVtxPos[1] = fESDVertex->GetY();
1213   fVtxPos[2] = fESDVertex->GetZ();
1214   Double_t CovMatx[6];
1215   fESDVertex->GetCovarianceMatrix(CovMatx);
1216   fVtxErr[0] = CovMatx[0];
1217   fVtxErr[1] = CovMatx[1];
1218   fVtxErr[2] = CovMatx[2];
1219   fVtxChi2 = fESDVertex->GetChi2();
1220   fVtxNDF = fESDVertex->GetNDF();
1221
1222   //Tracklets
1223   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1224
1225   //VZERO, ZDC
1226   AliESDVZERO *fV0data = esd->GetVZEROData();
1227   AliESDZDC *fZDCdata = esd->GetESDZDC();
1228
1229   fV0Adecision = fV0data->GetV0ADecision();
1230   fV0Cdecision = fV0data->GetV0CDecision();
1231   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1232   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1233
1234   fNLooseTracks = 0;
1235
1236   //Track loop - loose cuts
1237   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1238     AliESDtrack *trk = esd->GetTrack(itr);
1239     if( !trk ) continue;
1240
1241       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1242       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1243       if(trk->GetTPCNcls() < 20)continue;
1244       fNLooseTracks++;
1245   }//Track loop -loose cuts
1246
1247   Int_t nGoodTracks=0;
1248   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1249
1250   //Two Track loop
1251   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1252     AliESDtrack *trk = esd->GetTrack(itr);
1253     if( !trk ) continue;
1254
1255       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1256       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1257       if(trk->GetTPCNcls() < 70)continue;
1258       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1259       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1260       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1261       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1262       trk->GetImpactParameters(dca[0],dca[1]);
1263       if(TMath::Abs(dca[1]) > 2) continue;
1264       if(TMath::Abs(dca[1]) > 0.2) continue;
1265
1266       TrackIndex[nGoodTracks] = itr;
1267       nGoodTracks++;
1268       if(nGoodTracks > 2) break;
1269   }//Track loop
1270
1271   fJPsiESDTracks->Clear("C");
1272   if(nGoodTracks == 2){
1273           for(Int_t i=0; i<2; i++){
1274                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1275
1276                 AliExternalTrackParam cParam;
1277                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1278
1279                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
1280
1281                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1282                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1283                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1284
1285                 Double_t pos[3]={0,0,0};
1286                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1287                 else {
1288                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1289                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1290                      }
1291                 }
1292   if(!isMC) fJPsiTree ->Fill();
1293   }
1294
1295   nGoodTracks = 0;
1296   //Four track loop
1297   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1298     AliESDtrack *trk = esd->GetTrack(itr);
1299     if( !trk ) continue;
1300
1301       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1302       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1303       if(trk->GetTPCNcls() < 50)continue;
1304       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1305       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1306       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1307       trk->GetImpactParameters(dca[0],dca[1]);
1308       if(TMath::Abs(dca[1]) > 2) continue;
1309
1310       TrackIndex[nGoodTracks] = itr;
1311       nGoodTracks++;
1312       if(nGoodTracks > 4) break;
1313   }//Track loop
1314
1315   fPsi2sESDTracks->Clear("C");
1316   if(nGoodTracks == 4){
1317           for(Int_t i=0; i<4; i++){
1318                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1319
1320                 AliExternalTrackParam cParam;
1321                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1322
1323                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1324
1325                 fPIDMuon[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kMuon);
1326                 fPIDElectron[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kElectron);
1327                 fPIDPion[i] = fPIDResponse->NumberOfSigmasTPC(trk,AliPID::kPion);
1328
1329                 Double_t pos[3]={0,0,0};
1330                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1331                 else {
1332                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1333                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1334                      }
1335                 }
1336   if(!isMC) fPsi2sTree ->Fill();
1337   }
1338
1339   if(isMC){
1340         fJPsiTree ->Fill();
1341         fPsi2sTree ->Fill();
1342   }
1343
1344   PostData(1, fJPsiTree);
1345   PostData(2, fPsi2sTree);
1346
1347 }//RunESD
1348
1349
1350 //_____________________________________________________________________________
1351 void AliAnalysisTaskUpcPsi2s::RunESDMC(AliESDEvent* esd)
1352 {
1353
1354   AliTriggerAnalysis *fTrigAna = new AliTriggerAnalysis();
1355   fTrigAna->SetAnalyzeMC(isMC);
1356
1357   if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) > 1) fTriggerInputsMC[0] = kTRUE;
1358   if(fTrigAna->SPDFiredChips(esd,1,kFALSE,2) < 2) fTriggerInputsMC[0] = kFALSE;
1359
1360   fTriggerInputsMC[1] = esd->GetHeader()->IsTriggerInputFired("0OMU");
1361   fTriggerInputsMC[2] = esd->GetHeader()->IsTriggerInputFired("0VBA");
1362   fTriggerInputsMC[3] = esd->GetHeader()->IsTriggerInputFired("0VBC");
1363
1364   fGenPart->Clear("C");
1365
1366   AliMCEvent *mc = MCEvent();
1367   if(!mc) return;
1368
1369   Int_t nmc = 0;
1370   //loop over mc particles
1371   for(Int_t imc=0; imc<mc->GetNumberOfTracks(); imc++) {
1372     AliMCParticle *mcPart = (AliMCParticle*) mc->GetTrack(imc);
1373     if(!mcPart) continue;
1374
1375     if(mcPart->GetMother() >= 0) continue;
1376
1377     TParticle *part = (TParticle*) fGenPart->ConstructedAt(nmc++);
1378     part->SetMomentum(mcPart->Px(), mcPart->Py(), mcPart->Pz(), mcPart->E());
1379     part->SetPdgCode(mcPart->PdgCode());
1380     part->SetUniqueID(imc);
1381   }//loop over mc particles
1382
1383 }//RunESDMC
1384
1385
1386
1387 //_____________________________________________________________________________
1388 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
1389 {
1390
1391   cout<<"Analysis complete."<<endl;
1392 }//Terminate
1393