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[13] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","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",13,0.5,13.5);
 273   for (Int_t i = 0; i<13; 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[14] = {"Analyzed","Triggered","Vertex cut","V0 decision","Neutron ZDC cut","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",14,0.5,14.5);
 294   for (Int_t i = 0; i<14; 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) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
 360
 361   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0) fHistCentralTriggersPerRun->Fill(fRunNum);
 362
 363   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
 364
 365 PostData(3, fListTrig);
 366
 367 }
 368 //_____________________________________________________________________________
 369 void AliAnalysisTaskUpcPsi2s::RunAODhist()
 370 {
 371
 372   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 373
 374   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 375   Double_t muonMass = partMuon->Mass();
 376
 377   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 378   Double_t electronMass = partElectron->Mass();
 379
 380   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 381   Double_t pionMass = partPion->Mass();
 382
 383   //input event
 384   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 385   if(!aod) return;
 386
 387   fHistNeventsJPsi->Fill(1);
 388   fHistNeventsPsi2s->Fill(1);
 389
 390   //Trigger
 391   TString trigger = aod->GetFiredTriggerClasses();
 392
 393   if( !trigger.Contains("CCUP") ) return;
 394
 395   fHistNeventsJPsi->Fill(2);
 396   fHistNeventsPsi2s->Fill(2);
 397
 398   //primary vertex
 399   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 400   fVtxContrib = fAODVertex->GetNContributors();
 401   if(fVtxContrib < 2) return;
 402
 403   fHistNeventsJPsi->Fill(3);
 404   fHistNeventsPsi2s->Fill(3);
 405
 406   //VZERO, ZDC
 407   AliAODVZERO *fV0data = aod ->GetVZEROData();
 408   AliAODZDC *fZDCdata = aod->GetZDCData();
 409
 410   fV0Adecision = fV0data->GetV0ADecision();
 411   fV0Cdecision = fV0data->GetV0CDecision();
 412   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
 413
 414   fHistNeventsJPsi->Fill(4);
 415   fHistNeventsPsi2s->Fill(4);
 416
 417   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 418   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 419
 420   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
 421
 422   fHistNeventsJPsi->Fill(5);
 423   fHistNeventsPsi2s->Fill(5);
 424
 425   //Two tracks loop
 426   Int_t nGoodTracks = 0;
 427   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 428
 429   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
 430   Short_t qLepton[4], qPion[4];
 431   UInt_t nLepton=0, nPion=0, nHighPt=0;
 432   Double_t fRecTPCsignal[5];
 433   Int_t mass[3]={-1,-1,-1};
 434
 435
 436   //Four track loop
 437   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 438     AliAODTrack *trk = aod->GetTrack(itr);
 439     if( !trk ) continue;
 440     if(!(trk->TestFilterBit(1<<0))) continue;
 441
 442       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 443       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 444       if(trk->GetTPCNcls() < 50)continue;
 445       if(trk->Chi2perNDF() > 4)continue;
 446       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 447       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 448       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 449       delete trk_clone;
 450
 451       if(TMath::Abs(dca[1]) > 2) continue;
 452
 453       TrackIndex[nGoodTracks] = itr;
 454       nGoodTracks++;
 455
 456       if(nGoodTracks > 4) break;
 457   }//Track loop
 458
 459   nLepton=0; nPion=0; nHighPt=0;
 460   mass[0]= -1; mass[1]= -1, mass[2]= -1;
 461
 462   if(nGoodTracks == 4){
 463           fHistNeventsPsi2s->Fill(6);
 464           for(Int_t i=0; i<4; i++){
 465                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 466
 467                 if(trk->Pt() > 1){
 468                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 469                         qLepton[nLepton] = trk->Charge();
 470                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 471                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 472                                         mass[nLepton] = 0;
 473                                         }
 474                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 475                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 476                                         mass[nLepton] = 1;
 477                                         }
 478                         nLepton++;
 479                         }
 480                 else{
 481                         qPion[nPion] = trk->Charge();
 482                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
 483                         nPion++;
 484                         }
 485
 486                 if(nLepton > 2 || nPion > 2) break;
 487                 }
 488         if((nLepton == 2) && (nPion == 2)){
 489                 fHistNeventsPsi2s->Fill(7);
 490                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(8);
 491                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(9);
 492                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(10);
 493                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
 494                         fHistNeventsPsi2s->Fill(11);
 495                         if(mass[0] == mass[1]) {
 496                                 fHistNeventsPsi2s->Fill(12);
 497                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
 498                                 vDilepton = vLepton[0]+vLepton[1];
 499                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
 500                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
 501                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(13);
 502                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(14);
 503                                 }
 504                         }
 505                 }
 506   }
 507
 508   nGoodTracks = 0;
 509   //Two track loop
 510   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 511     AliAODTrack *trk = aod->GetTrack(itr);
 512     if( !trk ) continue;
 513     if(!(trk->TestFilterBit(1<<0))) 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(6);
 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(7);
 555                 if(qLepton[0]*qLepton[1] < 0){
 556                         fHistNeventsJPsi->Fill(8);
 557                         if(nHighPt > 0){
 558                                 fHistNeventsJPsi->Fill(9);
 559                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
 560                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(10);
 561                                 if(mass[0] == mass[1] && mass[0] != -1) {
 562                                         fHistNeventsJPsi->Fill(11);
 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(12);
 568                                                 }
 569                                         if(mass[0] == 1) {
 570                                                 fHistDiElectronMass->Fill(vCandidate.M());
 571                                                 fHistNeventsJPsi->Fill(13);
 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     if(!(trk->TestFilterBit(1<<0))) continue;
 653
 654       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 655       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 656       if(trk->GetTPCNcls() < 20)continue;
 657       fNLooseTracks++;
 658   }//Track loop -loose cuts
 659
 660   Int_t nGoodTracks=0;
 661   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 662
 663   //Two track loop
 664   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 665     AliAODTrack *trk = aod->GetTrack(itr);
 666     if( !trk ) continue;
 667     if(!(trk->TestFilterBit(1<<0))) continue;
 668
 669       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 670       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 671       if(trk->GetTPCNcls() < 70)continue;
 672       if(trk->Chi2perNDF() > 4)continue;
 673       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 674       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 675       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 676       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 677       delete trk_clone;
 678       if(TMath::Abs(dca[1]) > 2) continue;
 679       if(TMath::Abs(dca[0]) > 0.2) continue;
 680
 681       TrackIndex[nGoodTracks] = itr;
 682       nGoodTracks++;
 683
 684       if(nGoodTracks > 2) break;
 685   }//Track loop
 686
 687   if(nGoodTracks == 2){
 688
 689           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 690           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 691           Double_t muonMass = partMuon->Mass();
 692           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 693           Double_t electronMass = partElectron->Mass();
 694           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 695           Double_t pionMass = partPion->Mass();
 696
 697           Double_t KFcov[21];
 698           Double_t KFpar[6];
 699           Double_t KFmass = pionMass;
 700           Double_t fRecTPCsignal;
 701           AliKFParticle *KFpart[2];
 702           AliKFVertex *KFvtx = new AliKFVertex();
 703           KFvtx->SetField(aod->GetMagneticField());
 704
 705           for(Int_t i=0; i<2; i++){
 706                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 707
 708                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 709                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 710                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 711                 delete trk_clone;
 712
 713                 new((*fJPsiAODTracks)[i]) AliAODTrack(*trk);
 714                 ((AliAODTrack*)((*fJPsiAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
 715
 716                 trk->GetPosition(KFpar);
 717                 trk->PxPyPz(KFpar+3);
 718                 trk->GetCovMatrix(KFcov);
 719
 720                 if(trk->Pt() > 1){
 721                         fRecTPCsignal = trk->GetTPCsignal();
 722                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
 723                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
 724                         }
 725                 else KFmass = pionMass;
 726
 727                 KFpart[i] = new AliKFParticle();
 728                 KFpart[i]->SetField(aod->GetMagneticField());
 729                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
 730                 KFvtx->AddDaughter(*KFpart[i]);
 731
 732
 733                 Double_t pos[3]={0,0,0};
 734                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
 735                 parTrk->CopyFromVTrack((AliVTrack*) trk);
 736                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
 737                 else {
 738                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
 739                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
 740                      }
 741                 delete parTrk;
 742                 }
 743   fKfVtxPos[0]= KFvtx->GetX();
 744   fKfVtxPos[1]= KFvtx->GetY();
 745   fKfVtxPos[2]= KFvtx->GetZ();
 746   for(UInt_t i=0; i<2; i++)delete KFpart[i];
 747   delete KFvtx;
 748
 749   fJPsiTree ->Fill();
 750   }
 751
 752    nGoodTracks = 0;
 753    //Four track loop
 754   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 755     AliAODTrack *trk = aod->GetTrack(itr);
 756     if( !trk ) continue;
 757     if(!(trk->TestFilterBit(1<<0))) continue;
 758
 759       if(!(trk->GetStatus() & AliAODTrack::kTPCrefit) ) continue;
 760       if(!(trk->GetStatus() & AliAODTrack::kITSrefit) ) continue;
 761       if(trk->GetTPCNcls() < 50)continue;
 762       if(trk->Chi2perNDF() > 4)continue;
 763       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 764       AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 765       if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 766       delete trk_clone;
 767       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 768       if(TMath::Abs(dca[1]) > 2) continue;
 769
 770       TrackIndex[nGoodTracks] = itr;
 771       nGoodTracks++;
 772
 773       if(nGoodTracks > 4) break;
 774   }//Track loop
 775
 776
 777   if(nGoodTracks == 4){
 778
 779           TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 780           TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 781           Double_t muonMass = partMuon->Mass();
 782           TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 783           Double_t electronMass = partElectron->Mass();
 784           TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 785           Double_t pionMass = partPion->Mass();
 786
 787           Double_t KFcov[21];
 788           Double_t KFpar[6];
 789           Double_t KFmass = pionMass;
 790           Double_t fRecTPCsignal;
 791           AliKFParticle *KFpart[4];
 792           AliKFVertex *KFvtx = new AliKFVertex();
 793           KFvtx->SetField(aod->GetMagneticField());
 794
 795           for(Int_t i=0; i<4; i++){
 796                 AliAODTrack *trk = aod->GetTrack(TrackIndex[i]);
 797
 798                 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 799                 AliAODTrack* trk_clone=(AliAODTrack*)trk->Clone("trk_clone");
 800                 if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 801                 delete trk_clone;
 802
 803                 new((*fPsi2sAODTracks)[i]) AliAODTrack(*trk);
 804                 ((AliAODTrack*)((*fPsi2sAODTracks)[i]))->SetDCA(dca[0],dca[1]);//to get DCAxy trk->DCA(); to get DCAz trk->ZAtDCA();
 805
 806                 trk->GetPosition(KFpar);
 807                 trk->PxPyPz(KFpar+3);
 808                 trk->GetCovMatrix(KFcov);
 809
 810                 if(trk->Pt() > 1){
 811                         fRecTPCsignal = trk->GetTPCsignal();
 812                         if(fRecTPCsignal > 40 && fRecTPCsignal < 70) KFmass = muonMass;
 813                         if(fRecTPCsignal > 70 && fRecTPCsignal < 100)KFmass = electronMass;
 814                         }
 815                 else KFmass = pionMass;
 816
 817                 KFpart[i] = new AliKFParticle();
 818                 KFpart[i]->SetField(aod->GetMagneticField());
 819                 KFpart[i]->AliKFParticleBase::Initialize(KFpar,KFcov,(Int_t) trk->Charge(), KFmass);
 820                 KFvtx->AddDaughter(*KFpart[i]);
 821
 822                 Double_t pos[3]={0,0,0};
 823                 AliExternalTrackParam *parTrk = new AliExternalTrackParam();
 824                 parTrk->CopyFromVTrack((AliVTrack*) trk);
 825                 if(!parTrk->GetXYZAt(378,aod->GetMagneticField(),pos)) fTOFphi[i] = -666;
 826                 else {
 827                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
 828                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
 829                      }
 830                 delete parTrk;
 831                 }
 832   fKfVtxPos[0]= KFvtx->GetX();
 833   fKfVtxPos[1]= KFvtx->GetY();
 834   fKfVtxPos[2]= KFvtx->GetZ();
 835   for(UInt_t i=0; i<4; i++)delete KFpart[i];
 836   delete KFvtx;
 837   fPsi2sTree ->Fill();
 838   }
 839
 840   PostData(1, fJPsiTree);
 841   PostData(2, fPsi2sTree);
 842
 843 }//RunAOD
 844
 845 //_____________________________________________________________________________
 846 void AliAnalysisTaskUpcPsi2s::RunESDtrig()
 847 {
 848
 849   //input event
 850   AliESDEvent *esd = (AliESDEvent*) InputEvent();
 851   if(!esd) return;
 852
 853   fRunNum = esd ->GetRunNumber();
 854   //Trigger
 855   TString trigger = esd->GetFiredTriggerClasses();
 856
 857   if(trigger.Contains("CCUP4-B")) fHistUpcTriggersPerRun->Fill(fRunNum); //Upc triggers
 858
 859   if(trigger.Contains("CVLN_B2-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - synchronously downscaled
 860   if(trigger.Contains("CVLN_R1-B")) fHistCvlnTriggersPerRun->Fill(fRunNum); //CVLN triggers - randomly downscaled
 861
 862   if(esd->GetHeader()->IsTriggerInputFired("1ZED")) fHistZedTriggersPerRun->Fill(fRunNum); //1ZED trigger inputs
 863
 864    //MB, Central and SemiCentral triggers
 865   AliCentrality *centrality = esd->GetCentrality();
 866   UInt_t selectionMask = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 867
 868   //Double_t percentile = centrality->GetCentralityPercentile("V0M");
 869   Double_t percentile = centrality->GetCentralityPercentileUnchecked("V0M");
 870
 871   if(((selectionMask & AliVEvent::kMB) == AliVEvent::kMB) && percentile<=80 && percentile>=0) fHistMBTriggersPerRun->Fill(fRunNum);
 872
 873   if(((selectionMask & AliVEvent::kCentral) == AliVEvent::kCentral) && percentile<=6 && percentile>=0) fHistCentralTriggersPerRun->Fill(fRunNum);
 874
 875   if(((selectionMask & AliVEvent::kSemiCentral) == AliVEvent::kSemiCentral) && percentile<=50 && percentile>=15) fHistSemiCentralTriggersPerRun->Fill(fRunNum);
 876
 877
 878 PostData(3, fListTrig);
 879
 880 }
 881 //_____________________________________________________________________________
 882 void AliAnalysisTaskUpcPsi2s::RunESDhist()
 883 {
 884
 885   TDatabasePDG *pdgdat = TDatabasePDG::Instance();
 886
 887   TParticlePDG *partMuon = pdgdat->GetParticle( 13 );
 888   Double_t muonMass = partMuon->Mass();
 889
 890   TParticlePDG *partElectron = pdgdat->GetParticle( 11 );
 891   Double_t electronMass = partElectron->Mass();
 892
 893   TParticlePDG *partPion = pdgdat->GetParticle( 211 );
 894   Double_t pionMass = partPion->Mass();
 895
 896   //input event
 897   AliESDEvent *esd = (AliESDEvent*) InputEvent();
 898   if(!esd) return;
 899
 900   fHistNeventsJPsi->Fill(1);
 901   fHistNeventsPsi2s->Fill(1);
 902
 903   //Trigger
 904   TString trigger = esd->GetFiredTriggerClasses();
 905
 906   if( !trigger.Contains("CCUP") ) return;
 907
 908   fHistNeventsJPsi->Fill(2);
 909   fHistNeventsPsi2s->Fill(2);
 910
 911   //primary vertex
 912   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
 913   fVtxContrib = fESDVertex->GetNContributors();
 914   if(fVtxContrib < 2) return;
 915
 916   fHistNeventsJPsi->Fill(3);
 917   fHistNeventsPsi2s->Fill(3);
 918
 919   //VZERO, ZDC
 920   AliESDVZERO *fV0data = esd->GetVZEROData();
 921   AliESDZDC *fZDCdata = esd->GetESDZDC();
 922
 923   fV0Adecision = fV0data->GetV0ADecision();
 924   fV0Cdecision = fV0data->GetV0CDecision();
 925   if(fV0Adecision != AliESDVZERO::kV0Empty || fV0Cdecision != AliESDVZERO::kV0Empty) return;
 926
 927   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
 928   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
 929   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
 930
 931   fHistNeventsJPsi->Fill(4);
 932   fHistNeventsPsi2s->Fill(4);
 933
 934    Int_t nGoodTracks=0;
 935   //Two tracks loop
 936   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
 937
 938   TLorentzVector vLepton[4], vPion[4], vCandidate, vDilepton;
 939   Short_t qLepton[4], qPion[4];
 940   UInt_t nLepton=0, nPion=0, nHighPt=0;
 941   Double_t fRecTPCsignal[5];
 942   Int_t mass[3]={-1,-1,-1};
 943
 944  //Two Track loop
 945   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
 946     AliESDtrack *trk = esd->GetTrack(itr);
 947     if( !trk ) continue;
 948
 949       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 950       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 951       if(trk->GetTPCNcls() < 70)continue;
 952       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
 953       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
 954       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
 955       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
 956       trk->GetImpactParameters(dca[0],dca[1]);
 957       if(TMath::Abs(dca[1]) > 2) continue;
 958       if(TMath::Abs(dca[1]) > 0.2) continue;
 959
 960       TrackIndex[nGoodTracks] = itr;
 961       nGoodTracks++;
 962       if(nGoodTracks > 2) break;
 963   }//Track loop
 964
 965
 966   if(nGoodTracks == 2){
 967           fHistNeventsJPsi->Fill(5);
 968           for(Int_t i=0; i<2; i++){
 969                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
 970                 if(trk->Pt() > 1) nHighPt++;
 971                 fRecTPCsignal[nLepton] = trk->GetTPCsignal();
 972                 qLepton[nLepton] = trk->Charge();
 973                 if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
 974                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
 975                                 mass[nLepton] = 0;
 976                                 }
 977                 if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
 978                                 vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
 979                                 mass[nLepton] = 1;
 980                                 }
 981                 nLepton++;
 982                 }
 983         if(nLepton == 2){
 984                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsJPsi->Fill(6);
 985                 if(qLepton[0]*qLepton[1] < 0){
 986                         fHistNeventsJPsi->Fill(7);
 987                         if(nHighPt > 0){
 988                                 fHistNeventsJPsi->Fill(8);
 989                                 fHistTPCsignalJPsi->Fill(fRecTPCsignal[0],fRecTPCsignal[1]);
 990                                 if(nHighPt == 2) fHistNeventsJPsi->Fill(9);
 991                                 if(mass[0] == mass[1] && mass[0] != -1) {
 992                                         fHistNeventsJPsi->Fill(10);
 993                                         vCandidate = vLepton[0]+vLepton[1];
 994                                         if( vCandidate.M() > 2.8 && vCandidate.M() < 3.2) fHistDiLeptonPtJPsi->Fill(vLepton[0].Pt(),vLepton[1].Pt());
 995                                         if(mass[0] == 0) {
 996                                                 fHistDiMuonMass->Fill(vCandidate.M());
 997                                                 fHistNeventsJPsi->Fill(11);
 998                                                 }
 999                                         if(mass[0] == 1) {
1000                                                 fHistDiElectronMass->Fill(vCandidate.M());
1001                                                 fHistNeventsJPsi->Fill(12);
1002                                                 }
1003                                         }
1004                                 }
1005                         }
1006                 }
1007   }
1008   nGoodTracks = 0; nLepton=0; nPion=0; nHighPt=0;
1009   mass[0]= -1; mass[1]= -1, mass[2]= -1;
1010
1011     //Four Track loop
1012   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1013     AliESDtrack *trk = esd->GetTrack(itr);
1014     if( !trk ) continue;
1015
1016       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1017       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1018       if(trk->GetTPCNcls() < 50)continue;
1019       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1020       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1021       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1022       trk->GetImpactParameters(dca[0],dca[1]);
1023       if(TMath::Abs(dca[1]) > 2) continue;
1024
1025       TrackIndex[nGoodTracks] = itr;
1026       nGoodTracks++;
1027       if(nGoodTracks > 4) break;
1028   }//Track loop
1029
1030   if(nGoodTracks == 4){
1031           fHistNeventsPsi2s->Fill(5);
1032           for(Int_t i=0; i<4; i++){
1033                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1034
1035                 if(trk->Pt() > 1){
1036                         fRecTPCsignal[nLepton] = trk->GetTPCsignal();
1037                         qLepton[nLepton] = trk->Charge();
1038                         if(fRecTPCsignal[nLepton] > 40 && fRecTPCsignal[nLepton] < 70){
1039                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), muonMass);
1040                                         mass[nLepton] = 0;
1041                                         }
1042                         if(fRecTPCsignal[nLepton] > 70 && fRecTPCsignal[nLepton] < 100){
1043                                         vLepton[nLepton].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), electronMass);
1044                                         mass[nLepton] = 1;
1045                                         }
1046                         nLepton++;
1047                         }
1048                 else{
1049                         qPion[nPion] = trk->Charge();
1050                         vPion[nPion].SetPtEtaPhiM(trk->Pt(), trk->Eta(), trk->Phi(), pionMass);
1051                         nPion++;
1052                         }
1053
1054                 if(nLepton > 2 || nPion > 2) break;
1055                 }
1056         if((nLepton == 2) && (nPion == 2)){
1057                 fHistNeventsPsi2s->Fill(6);
1058                 if(qLepton[0]*qLepton[1] > 0) fHistNeventsPsi2s->Fill(7);
1059                 if(qPion[0]*qPion[1] > 0) fHistNeventsPsi2s->Fill(8);
1060                 if((qLepton[0]*qLepton[1] > 0) && (qPion[0]*qPion[1] > 0)) fHistNeventsPsi2s->Fill(9);
1061                 if((qLepton[0]*qLepton[1] < 0) && (qPion[0]*qPion[1] < 0)){
1062                         fHistNeventsPsi2s->Fill(10);
1063                         if(mass[0] == mass[1]) {
1064                                 fHistNeventsPsi2s->Fill(11);
1065                                 vCandidate = vLepton[0]+vLepton[1]+vPion[0]+vPion[1];
1066                                 vDilepton = vLepton[0]+vLepton[1];
1067                                 fHistPsi2sMassVsPt->Fill(vCandidate.M(),vCandidate.Pt());
1068                                 if(vCandidate.Pt() < 0.15) fHistPsi2sMassCoherent->Fill(vCandidate.M());
1069                                 if(mass[0] == 0) fHistNeventsPsi2s->Fill(12);
1070                                 if(mass[0] == 1) fHistNeventsPsi2s->Fill(13);
1071                                 }
1072                         }
1073                 }
1074   }
1075
1076   PostData(4, fListHist);
1077
1078 }
1079
1080 //_____________________________________________________________________________
1081 void AliAnalysisTaskUpcPsi2s::RunESDtree()
1082 {
1083
1084   //input event
1085   AliESDEvent *esd = (AliESDEvent*) InputEvent();
1086   if(!esd) return;
1087
1088   //input data
1089   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
1090   fDataFilnam->Clear();
1091   fDataFilnam->SetString(filnam);
1092   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
1093   fRunNum = esd->GetRunNumber();
1094
1095    //Trigger
1096   TString trigger = esd->GetFiredTriggerClasses();
1097
1098   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
1099   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
1100   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
1101
1102   Bool_t isTriggered = kFALSE;
1103   for(Int_t i=0; i<ntrg; i++) {
1104     if( fTrigger[i] ) isTriggered = kTRUE;
1105   }
1106   if( !isTriggered ) return;
1107
1108   //trigger inputs
1109   fL0inputs = esd->GetHeader()->GetL0TriggerInputs();
1110   fL1inputs = esd->GetHeader()->GetL1TriggerInputs();
1111
1112   //TOF trigger info (0OMU)
1113   fTOFtrig1 = esd->GetHeader()->IsTriggerInputFired("0OMU");
1114   fTOFtrig2 = esd->GetHeader()->GetActiveTriggerInputs().Contains("0OMU") ? ((esd->GetHeader()) ? esd->GetHeader()->IsTriggerInputFired("0OMU") : kFALSE) : TESTBIT(esd->GetHeader()->GetL0TriggerInputs(), (kFALSE) ? 21 : 9);
1115
1116   //Event identification
1117   fPerNum = esd->GetPeriodNumber();
1118   fOrbNum = esd->GetOrbitNumber();
1119   fBCrossNum = esd->GetBunchCrossNumber();
1120
1121   //primary vertex
1122   AliESDVertex *fESDVertex = (AliESDVertex*) esd->GetPrimaryVertex();
1123   fVtxContrib = fESDVertex->GetNContributors();
1124   fVtxPos[0] = fESDVertex->GetX();
1125   fVtxPos[1] = fESDVertex->GetY();
1126   fVtxPos[2] = fESDVertex->GetZ();
1127   Double_t CovMatx[6];
1128   fESDVertex->GetCovarianceMatrix(CovMatx);
1129   fVtxErr[0] = CovMatx[0];
1130   fVtxErr[1] = CovMatx[1];
1131   fVtxErr[2] = CovMatx[2];
1132   fVtxChi2 = fESDVertex->GetChi2();
1133   fVtxNDF = fESDVertex->GetNDF();
1134
1135   //Tracklets
1136   fNtracklets = esd->GetMultiplicity()->GetNumberOfTracklets();
1137
1138   //VZERO, ZDC
1139   AliESDVZERO *fV0data = esd->GetVZEROData();
1140   AliESDZDC *fZDCdata = esd->GetESDZDC();
1141
1142   fV0Adecision = fV0data->GetV0ADecision();
1143   fV0Cdecision = fV0data->GetV0CDecision();
1144   fZDCAenergy = fZDCdata->GetZN2TowerEnergy()[0];
1145   fZDCCenergy = fZDCdata->GetZN1TowerEnergy()[0];
1146
1147   fNLooseTracks = 0;
1148
1149   //Track loop - loose cuts
1150   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1151     AliESDtrack *trk = esd->GetTrack(itr);
1152     if( !trk ) continue;
1153
1154       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1155       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1156       if(trk->GetTPCNcls() < 20)continue;
1157       fNLooseTracks++;
1158   }//Track loop -loose cuts
1159
1160   Int_t nGoodTracks=0;
1161   Int_t TrackIndex[5] = {-1,-1,-1,-1,-1};
1162
1163   //Two Track loop
1164   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1165     AliESDtrack *trk = esd->GetTrack(itr);
1166     if( !trk ) continue;
1167
1168       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1169       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1170       if(trk->GetTPCNcls() < 70)continue;
1171       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1172       if((!trk->HasPointOnITSLayer(0))&&(!trk->HasPointOnITSLayer(1))) continue;
1173       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1174       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1175       trk->GetImpactParameters(dca[0],dca[1]);
1176       if(TMath::Abs(dca[1]) > 2) continue;
1177       if(TMath::Abs(dca[1]) > 0.2) continue;
1178
1179       TrackIndex[nGoodTracks] = itr;
1180       nGoodTracks++;
1181       if(nGoodTracks > 2) break;
1182   }//Track loop
1183
1184   if(nGoodTracks == 2){
1185           for(Int_t i=0; i<2; i++){
1186                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1187
1188                 AliExternalTrackParam cParam;
1189                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1190
1191                 new((*fJPsiESDTracks)[i]) AliESDtrack(*trk);
1192
1193                 Double_t pos[3]={0,0,0};
1194                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1195                 else {
1196                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1197                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1198                      }
1199                 }
1200   fJPsiTree ->Fill();
1201   }
1202
1203   nGoodTracks = 0;
1204   //Four track loop
1205   for(Int_t itr=0; itr<esd ->GetNumberOfTracks(); itr++) {
1206     AliESDtrack *trk = esd->GetTrack(itr);
1207     if( !trk ) continue;
1208
1209       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
1210       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
1211       if(trk->GetTPCNcls() < 50)continue;
1212       if(trk->GetTPCchi2()/trk->GetTPCNcls() > 4)continue;
1213       Float_t dca[2] = {0.0,0.0}; AliExternalTrackParam cParam;
1214       if(!trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam)) continue;
1215       trk->GetImpactParameters(dca[0],dca[1]);
1216       if(TMath::Abs(dca[1]) > 2) continue;
1217
1218       TrackIndex[nGoodTracks] = itr;
1219       nGoodTracks++;
1220       if(nGoodTracks > 4) break;
1221   }//Track loop
1222
1223   if(nGoodTracks == 4){
1224           for(Int_t i=0; i<4; i++){
1225                 AliESDtrack *trk = esd->GetTrack(TrackIndex[i]);
1226
1227                 AliExternalTrackParam cParam;
1228                 trk->RelateToVertex(fESDVertex, esd->GetMagneticField(),300.,&cParam);// to get trk->GetImpactParameters(DCAxy,DCAz);
1229
1230                 new((*fPsi2sESDTracks)[i]) AliESDtrack(*trk);
1231
1232                 Double_t pos[3]={0,0,0};
1233                 if(!trk->GetXYZAt(378,esd->GetMagneticField(),pos)) fTOFphi[i] = -666;
1234                 else {
1235                      fTOFphi[i] =  TMath::ATan2(pos[1],pos[0])*TMath::RadToDeg();
1236                      if(fTOFphi[i] < 0) fTOFphi[i]+=(2*TMath::Pi()*TMath::RadToDeg());
1237                      }
1238                 }
1239   fPsi2sTree ->Fill();
1240   }
1241
1242   PostData(1, fJPsiTree);
1243   PostData(2, fPsi2sTree);
1244
1245 }//RunESD
1246
1247 //_____________________________________________________________________________
1248 void AliAnalysisTaskUpcPsi2s::Terminate(Option_t *)
1249 {
1250
1251   cout<<"Analysis complete."<<endl;
1252 }//Terminate
1253