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