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