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